Nothing
R/tabmulti.R
In R4VN: Health Data Analysis and Publication-Ready Reporting
Defines functions
print.r4vn_tabmulti
tabmulti
Documented in
print.r4vn_tabmulti
tabmulti
# tabmulti_complete.R
# Rebuilt from the user's tabmulti(1).R and aligned with tab(10).R API/rendering.
# LASSO fix: lasso_lambda is validated explicitly and never passed to match.arg().
#' Compare Multivariable Model-Building Strategies
#'
#' Builds and compares variable-selection strategies for a binary outcome.
#' Each strategy selects complete variables or terms, after which the selected
#' model is refitted using ordinary logistic regression or modified Poisson
#' regression so that conventional OR, RR, or PR estimates, 95% confidence
#' intervals, and p-values can be reported.
#'
#' @usage
#' tabmulti(data = NULL, vars = NULL, by = NULL,
#' methods = c("full", "forward", "backward", "purposeful"),
#' digit = 1, p_digit = 3, effect_digit = 2, global = FALSE,
#' pvalue = TRUE, rvrow = NULL, bold_p = TRUE, p_bold = 0.05,
#' or = FALSE, rr = FALSE, pr = FALSE, event = NULL,
#' criterion = c("AIC", "BIC"), force = NULL, entry = 0.20,
#' stay = 0.05, confounding = 0.10, lasso_lambda = "lambda.1se",
#' bma_pip = 0.50, max_subset_vars = 15L, max_subset_models = 100000L,
#' template = c("journal", "clean", "minimal"), append = NULL,
#' file = NULL, raw = FALSE, name = FALSE, title = NULL, show = TRUE)
#'
#' @param data Optional data frame. When omitted or \code{NULL}, the active
#' data frame set by \code{usedf()} or \code{opendata(..., active = TRUE)}
#' is used.
#' @param vars A variable specification created by \code{vars()}.
#' @param by Binary outcome supplied without quotation marks.
#' @param methods Model strategies: \code{"full"}, \code{"forward"},
#' \code{"backward"}, \code{"purposeful"}, \code{"lasso"},
#' \code{"bma"}, or \code{"best"}.
#' @param digit Retained for API consistency with \code{tab()}.
#' @param p_digit Number of decimal places for p-values.
#' @param effect_digit Number of decimal places for effect estimates and
#' confidence limits.
#' @param global Logical. Display global likelihood-ratio p-values.
#' @param pvalue Logical. Display coefficient p-value columns.
#' @param rvrow Categorical variables whose displayed level order should be
#' reversed. This does not change model reference categories.
#' @param bold_p Logical. Bold p-values smaller than \code{p_bold}.
#' @param p_bold Threshold used when \code{bold_p = TRUE}.
#' @param or Logical. Report odds ratios from logistic regression.
#' @param rr Logical. Report risk ratios from modified Poisson regression.
#' @param pr Logical. Report prevalence ratios from modified Poisson regression.
#' Exactly one of \code{or}, \code{rr}, and \code{pr} must be
#' \code{TRUE}.
#' @param event Event level. The last observed outcome level is used when omitted.
#' @param criterion Selection criterion, \code{"AIC"} or \code{"BIC"}.
#' @param force Variables forced into every selected model. Accepts
#' \code{vars(...)}, \code{c(...)}, a character vector, \code{TRUE}, or
#' \code{"ALL"}.
#' @param entry Univariate entry threshold for purposeful selection.
#' @param stay Multivariable retention threshold for purposeful selection.
#' @param confounding Relative coefficient-change threshold for identifying a
#' confounder during purposeful selection.
#' @param lasso_lambda Either \code{"lambda.min"} or \code{"lambda.1se"}.
#' @param bma_pip Posterior inclusion-probability threshold used by the
#' BIC-weighted BMA strategy.
#' @param max_subset_vars Maximum number of candidate variables for exhaustive
#' subset methods.
#' @param max_subset_models Maximum number of subset models to evaluate.
#' @param template HTML style: \code{"journal"}, \code{"clean"}, or
#' \code{"minimal"}.
#' @param append Optional previous \code{r4vn_tabmulti} object or existing HTML path.
#' @param file Optional output HTML path.
#' @param raw Logical. Retain unformatted coefficients and selection details.
#' @param name Logical. Display original variable names beside labels.
#' @param title Optional table title.
#' @param show Logical. Open the HTML table in the Viewer or browser.
#'
#' @details
#' Available strategies are:
#' \itemize{
#' \item \code{full}: include every candidate variable;
#' \item \code{forward}: forward stepwise selection using AIC or BIC;
#' \item \code{backward}: backward stepwise selection from the full model;
#' \item \code{purposeful}: univariate screening followed by significance
#' and confounding assessment;
#' \item \code{lasso}: selection with \code{glmnet::cv.glmnet()}, followed
#' by ordinary-model refitting; requires the suggested package \pkg{glmnet};
#' \item \code{bma}: BIC-weighted subset averaging and inclusion-probability
#' thresholding;
#' \item \code{best}: select the subset with the smallest AIC or BIC.
#' }
#'
#' All strategies use the same complete-case sample. Diagnostic rows include
#' sample size, events, number of variables and parameters, AIC, BIC,
#' pseudo-R-squared measures, goodness-of-fit tests, AUC where applicable, and
#' the coefficient-level VIF range. Exhaustive methods grow exponentially with
#' the number of candidate variables.
#'
#' @return Invisibly returns an object of class \code{r4vn_tabmulti}. Important
#' components include \code{data}, \code{selected}, \code{models},
#' \code{diagnostics}, \code{file}, \code{html}, and
#' \code{table_html}.
#'
#' @seealso \code{\link{vars}}, \code{\link{tab}}, and
#' \code{\link{tabexport}}.
#' @family R4VN tables
#'
#' @examples
#' set.seed(2026)
#' n <- 180
#' dat <- data.frame(
#' age = round(rnorm(n, 45, 12)),
#' sex = factor(sample(c("Female", "Male"), n, TRUE)),
#' bmi = round(rnorm(n, 23, 3), 1),
#' smoking = factor(sample(c("No", "Yes"), n, TRUE,
#' prob = c(0.70, 0.30))),
#' education = factor(sample(c("Primary", "Secondary", "College"),
#' n, TRUE))
#' )
#' lp <- -3.1 + 0.045 * dat$age + 0.11 * (dat$bmi - 23) +
#' 0.45 * (dat$sex == "Male") + 0.70 * (dat$smoking == "Yes")
#' dat$hypertension <- factor(
#' rbinom(n, 1, plogis(lp)),
#' levels = c(0, 1), labels = c("No", "Yes")
#' )
#'
#' models <- tabmulti(
#' dat,
#' vars = vars(c.age, b2.sex, c.bmi, b2.smoking, b2.education),
#' by = hypertension,
#' methods = c("full", "backward"),
#' or = TRUE,
#' event = "Yes",
#' criterion = "AIC",
#' global = TRUE,
#' show = FALSE
#' )
#' models$selected
#' models$diagnostics
#'
#' \donttest{
#' # LASSO requires the suggested package glmnet.
#' if (requireNamespace("glmnet", quietly = TRUE)) {
#' models_lasso <- tabmulti(
#' dat,
#' vars = vars(c.age, b2.sex, c.bmi, b2.smoking, b2.education),
#' by = hypertension,
#' methods = c("full", "lasso"),
#' or = TRUE,
#' event = "Yes",
#' show = FALSE
#' )
#' }
#' }
#'
#' # Extended usage examples
#' \donttest{
#' set.seed(2026)
#' n <- 400
#' d <- data.frame(
#' sex = factor(sample(c("Female", "Male"), n, TRUE)),
#' age = rnorm(n, 45, 12),
#' bmi = rnorm(n, 23, 3),
#' smoking = factor(sample(c("No", "Yes"), n, TRUE)),
#' education = factor(sample(c("Primary", "Secondary", "College"), n, TRUE))
#' )
#' lp <- -3.2 + 0.045 * d$age + 0.10 * (d$bmi - 23) +
#' 0.45 * (d$sex == "Male") + 0.70 * (d$smoking == "Yes")
#' d$outcome <- factor(rbinom(n, 1, plogis(lp)),
#' levels = 0:1, labels = c("No", "Yes"))
#'
#' # Full multivariable logistic model
#' m1 <- tabmulti(
#' d,
#' vars = vars(b2.sex, c.age, c.bmi, b2.smoking, b2.education),
#' by = outcome,
#' methods = "full",
#' or = TRUE,
#' event = "Yes",
#' show = FALSE
#' )
#'
#' # Compare several model-building strategies in one table
#' m2 <- tabmulti(
#' d,
#' vars = vars(b2.sex, c.age, c.bmi, b2.smoking, b2.education),
#' by = outcome,
#' methods = c("full", "forward", "backward", "purposeful"),
#' or = TRUE,
#' event = "Yes",
#' criterion = "AIC",
#' global = TRUE,
#' show = FALSE
#' )
#' m2$selected
#' m2$diagnostics
#'
#' # Force variables into every selected model and tune purposeful selection
#' tabmulti(
#' d,
#' vars = vars(b2.sex, c.age, c.bmi, b2.smoking, b2.education),
#' by = outcome,
#' methods = "purposeful",
#' force = vars(age, sex),
#' entry = 0.20, stay = 0.05, confounding = 0.10,
#' or = TRUE, event = "Yes", show = FALSE
#' )
#'
#' # Modified Poisson models for RR or PR
#' tabmulti(d, vars = vars(b2.sex, c.age, c.bmi, b2.smoking),
#' by = outcome, methods = "full", rr = TRUE,
#' event = "Yes", show = FALSE)
#' tabmulti(d, vars = vars(b2.sex, c.age, c.bmi, b2.smoking),
#' by = outcome, methods = "full", pr = TRUE,
#' event = "Yes", show = FALSE)
#'
#' # Display and output controls
#' tabmulti(d, vars = vars(b2.sex, c.age, c.bmi, b2.smoking),
#' by = outcome, methods = c("full", "backward"),
#' or = TRUE, event = "Yes", rvrow = vars(smoking),
#' pvalue = TRUE, bold_p = TRUE, p_bold = 0.05,
#' template = "clean", raw = TRUE, name = TRUE,
#' title = "Model-building comparison", show = FALSE)
#'
#' # Active-data syntax
#' usedf(d)
#' tabmulti(vars = vars(b2.sex, c.age, c.bmi, b2.smoking),
#' by = outcome, methods = "full", or = TRUE,
#' event = "Yes", show = FALSE)
#'
#' # LASSO requires glmnet; BMA/best are exhaustive and suit fewer candidates
#' if (requireNamespace("glmnet", quietly = TRUE)) {
#' tabmulti(d, vars = vars(b2.sex, c.age, c.bmi, b2.smoking),
#' by = outcome, methods = c("full", "lasso"), or = TRUE,
#' event = "Yes", lasso_lambda = "lambda.1se", show = FALSE)
#' }
#' tabmulti(d, vars = vars(b2.sex, c.age, c.bmi, b2.smoking),
#' by = outcome, methods = c("best", "bma"), or = TRUE,
#' event = "Yes", criterion = "BIC", bma_pip = 0.50,
#' max_subset_vars = 10, show = FALSE)
#' }
#' @export
tabmulti <- function(data = NULL, vars = NULL, by = NULL,
methods = c("full", "forward", "backward", "purposeful"),
digit = 1, p_digit = 3, effect_digit = 2,
global = FALSE, pvalue = TRUE, rvrow = NULL,
bold_p = TRUE, p_bold = 0.05,
or = FALSE, rr = FALSE, pr = FALSE, event = NULL,
criterion = c("AIC", "BIC"), force = NULL,
entry = 0.20, stay = 0.05, confounding = 0.10,
lasso_lambda = "lambda.1se", bma_pip = 0.50,
max_subset_vars = 15L, max_subset_models = 100000L,
template = c("journal", "clean", "minimal"),
append = NULL, file = NULL, raw = FALSE,
name = FALSE, title = NULL, show = TRUE) {
# Use active data when `data` is omitted.
if (is.null(data)) data <- .r4vn_get_active()
# ----------------------------- validation -----------------------------
if (!is.data.frame(data)) stop("`data` must be a data frame.", call. = FALSE)
if (!inherits(vars, "r4vn_vars")) stop("`vars` must be created using `vars()`.", call. = FALSE)
vars <- .r4vn_resolve_vars_input(
vars, data = data, arg = "vars", default_type = "auto", strict = TRUE
)
validate_integer <- function(x, arg) {
if (!is.numeric(x) || length(x) != 1L || is.na(x) || x < 0 || x != floor(x))
stop(sprintf("`%s` must be a single non-negative integer.", arg), call. = FALSE)
}
validate_flag <- function(x, arg) {
if (!is.logical(x) || length(x) != 1L || is.na(x))
stop(sprintf("`%s` must be TRUE or FALSE.", arg), call. = FALSE)
}
validate_probability <- function(x, arg) {
if (!is.numeric(x) || length(x) != 1L || is.na(x) || x < 0 || x > 1)
stop(sprintf("`%s` must be between 0 and 1.", arg), call. = FALSE)
}
validate_integer(digit, "digit")
validate_integer(p_digit, "p_digit")
validate_integer(effect_digit, "effect_digit")
validate_integer(max_subset_vars, "max_subset_vars")
validate_integer(max_subset_models, "max_subset_models")
for (a in c("global", "pvalue", "bold_p", "or", "rr", "pr", "raw", "name", "show"))
validate_flag(get(a), a)
validate_probability(p_bold, "p_bold")
validate_probability(entry, "entry")
validate_probability(stay, "stay")
validate_probability(confounding, "confounding")
validate_probability(bma_pip, "bma_pip")
template <- match.arg(template)
criterion <- match.arg(criterion)
allowed_methods <- c("full", "forward", "backward", "purposeful", "lasso", "bma", "best")
methods <- unique(tolower(as.character(methods)))
invalid <- setdiff(methods, allowed_methods)
if (!length(methods)) stop("`methods` must contain at least one method.", call. = FALSE)
if (length(invalid)) stop(sprintf("Unsupported methods: %s.", paste(invalid, collapse = ", ")), call. = FALSE)
# Do not use match.arg() here: exact validation avoids the previous ambiguous failure.
lasso_lambda <- as.character(lasso_lambda)[1L]
if (is.na(lasso_lambda) || !lasso_lambda %in% c("lambda.min", "lambda.1se"))
stop("`lasso_lambda` must be exactly 'lambda.min' or 'lambda.1se'.", call. = FALSE)
effects <- c(OR = or, RR = rr, PR = pr)
if (sum(effects) != 1L) stop("Exactly one of `or`, `rr`, or `pr` must be TRUE.", call. = FALSE)
effect_type <- names(effects)[effects][1L]
by_expr <- substitute(by)
if (!is.symbol(by_expr)) stop("`by` must be a single variable name.", call. = FALSE)
by_name <- as.character(by_expr)
duplicated_variables <- unique(vars$variable[duplicated(vars$variable)])
if (length(duplicated_variables))
stop(sprintf("Duplicated variable specification: %s.", paste(duplicated_variables, collapse = ", ")), call. = FALSE)
absent <- setdiff(c(vars$variable, by_name), names(data))
if (length(absent)) stop(sprintf("Variables not found in `data`: %s.", paste(absent, collapse = ", ")), call. = FALSE)
if (by_name %in% vars$variable) stop("The outcome variable cannot also appear in `vars`.", call. = FALSE)
# ----------------------------- helpers -----------------------------
escape_html <- 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)
}
get_label <- function(x, variable) {
z <- attr(x, "label", exact = TRUE)
if (is.null(z) || !length(z) || is.na(z[1L]) || !nzchar(as.character(z[1L]))) variable else as.character(z[1L])
}
display_label <- function(label, variable) {
if (!isTRUE(name)) return(escape_html(label))
paste0(escape_html(label), " <span class=\"variable-code\">[", escape_html(variable), "]</span>")
}
get_levels <- function(x) {
observed <- x[!is.na(x)]
values <- if (is.factor(x)) levels(x) else if (is.logical(x)) c(FALSE, TRUE) else {
z <- unique(observed)
if (is.numeric(z)) sort(z) else z
}
values[values %in% observed]
}
format_number <- function(x, digits = effect_digit) {
if (!length(x) || is.na(x) || !is.finite(x)) return("")
formatC(x, format = "f", digits = digits, big.mark = ",")
}
format_p_plain <- function(x) {
if (!length(x) || is.na(x) || !is.finite(x)) return("")
limit <- 10^(-p_digit)
if (x < limit) paste0("<", formatC(limit, format = "f", digits = p_digit))
else formatC(x, format = "f", digits = p_digit)
}
format_p_html <- function(x) {
z <- escape_html(format_p_plain(x))
if (nzchar(z) && isTRUE(bold_p) && is.finite(x) && x < p_bold) paste0("<strong>", z, "</strong>") else z
}
format_effect <- function(est, low, high) {
if (!all(is.finite(c(est, low, high)))) return("")
paste0(format_number(est), " (", format_number(low), " - ", format_number(high), ")")
}
make_formula <- function(variables) stats::reformulate(variables, response = ".outcome")
robust_vcov_poisson <- function(fit) {
X <- stats::model.matrix(fit)
mu <- stats::fitted(fit)
score <- fit$y - mu
bread <- tryCatch(solve(crossprod(X, X * as.vector(mu))), error = function(e) NULL)
if (is.null(bread)) return(NULL)
meat <- crossprod(X, X * as.vector(score^2))
out <- bread %*% meat %*% bread
dimnames(out) <- list(colnames(X), colnames(X))
out
}
# Parse force with the same accepted styles as tab(adjusted/multi).
parse_force <- function(expression) {
all_names <- vars$variable
if (identical(expression, quote(NULL))) return(character())
if (is.logical(expression) && length(expression) == 1L)
return(if (isTRUE(expression)) all_names else character())
if (is.character(expression)) {
if (length(expression) == 1L && toupper(expression) == "ALL") return(all_names)
return(expression)
}
if (is.symbol(expression)) {
object_name <- as.character(expression)
if (toupper(object_name) == "ALL") return(all_names)
value <- tryCatch(get(object_name, envir = parent.frame(2L)), error = function(e) NULL)
if (inherits(value, "r4vn_vars")) return(value$variable)
if (is.character(value)) return(if (length(value) == 1L && toupper(value) == "ALL") all_names else value)
if (is.logical(value) && length(value) == 1L && isTRUE(value)) return(all_names)
return(object_name)
}
if (is.call(expression) && as.character(expression[[1L]]) == "vars") {
value <- eval(expression, envir = parent.frame(2L))
if (!inherits(value, "r4vn_vars")) stop("`force = vars(...)` is invalid.", call. = FALSE)
return(value$variable)
}
if (is.call(expression) && as.character(expression[[1L]]) == "c") {
items <- as.list(expression)[-1L]
if (all(vapply(items, is.symbol, logical(1)))) return(vapply(items, as.character, character(1)))
value <- eval(expression, envir = parent.frame(2L))
if (is.character(value)) return(value)
}
value <- tryCatch(eval(expression, envir = parent.frame(2L)), error = function(e) NULL)
if (inherits(value, "r4vn_vars")) return(value$variable)
if (is.character(value)) return(value)
stop("`force` must be NULL, TRUE, 'ALL', vars(...), c(...), or a character vector.", call. = FALSE)
}
forced <- unique(parse_force(substitute(force)))
invalid_forced <- setdiff(forced, vars$variable)
if (length(invalid_forced)) stop(sprintf("Forced variables not listed in `vars`: %s.", paste(invalid_forced, collapse = ", ")), call. = FALSE)
# Parse rvrow using the same conventions as tab(). rvrow changes display
# order only; it never changes the reference category used in a model.
parse_rvrow <- function(expression) {
categorical <- vars$variable[vars$type == "categorical"]
if (identical(expression, quote(NULL))) return(character())
if (is.logical(expression) && length(expression) == 1L)
return(if (isTRUE(expression)) categorical else character())
if (is.character(expression)) return(expression)
if (is.symbol(expression)) {
object_name <- as.character(expression)
value <- tryCatch(get(object_name, envir = parent.frame(2L)), error = function(e) NULL)
if (inherits(value, "r4vn_vars")) return(value$variable)
if (is.character(value)) return(value)
if (is.logical(value) && length(value) == 1L) return(if (isTRUE(value)) categorical else character())
return(object_name)
}
if (is.call(expression) && as.character(expression[[1L]]) %in% c("vars", "c")) {
if (as.character(expression[[1L]]) == "vars") {
value <- eval(expression, envir = parent.frame(2L))
if (inherits(value, "r4vn_vars")) return(value$variable)
} else {
items <- as.list(expression)[-1L]
if (all(vapply(items, is.symbol, logical(1)))) return(vapply(items, as.character, character(1)))
value <- tryCatch(eval(expression, envir = parent.frame(2L)), error = function(e) NULL)
if (is.character(value)) return(value)
}
}
value <- tryCatch(eval(expression, envir = parent.frame(2L)), error = function(e) NULL)
if (inherits(value, "r4vn_vars")) return(value$variable)
if (is.character(value)) return(value)
stop("`rvrow` must be NULL, TRUE, vars(...), c(...), a variable name, or a character vector.", call. = FALSE)
}
reverse_rows <- unique(parse_rvrow(substitute(rvrow)))
invalid_reverse <- setdiff(reverse_rows, vars$variable[vars$type == "categorical"])
if (length(invalid_reverse)) stop(sprintf("Invalid categorical variables in `rvrow`: %s.", paste(invalid_reverse, collapse = ", ")), call. = FALSE)
# ----------------------------- outcome/data -----------------------------
y <- data[[by_name]]
original_y_levels <- get_levels(y)
if (length(original_y_levels) != 2L) stop("`by` must have exactly two observed levels.", call. = FALSE)
event_level <- if (is.null(event)) utils::tail(original_y_levels, 1L) else as.character(event)[1L]
if (!event_level %in% as.character(original_y_levels)) stop("`event` is not an observed outcome level.", call. = FALSE)
model_data <- data[c(by_name, vars$variable)]
model_data <- model_data[stats::complete.cases(model_data), , drop = FALSE]
if (!nrow(model_data)) stop("No complete observations remain for the requested variables.", call. = FALSE)
model_data$.outcome <- as.integer(as.character(model_data[[by_name]]) == event_level)
if (length(unique(model_data$.outcome)) != 2L) stop("The complete-case data do not contain both outcome groups.", call. = FALSE)
# Preserve metadata and references exactly as declared by vars().
meta <- vars
for (i in seq_len(nrow(meta))) {
v <- meta$variable[i]
if (meta$type[i] == "categorical") {
lev <- get_levels(data[[v]])
ref_index <- meta$reference_index[i]
if (is.na(ref_index) || ref_index < 1L || ref_index > length(lev))
stop(sprintf("Invalid reference index for `%s`.", v), call. = FALSE)
ref <- as.character(lev[ref_index])
model_data[[v]] <- factor(model_data[[v]], levels = lev)
model_data[[v]] <- stats::relevel(model_data[[v]], ref = ref)
} else {
model_data[[v]] <- suppressWarnings(as.numeric(model_data[[v]]))
if (anyNA(model_data[[v]])) stop(sprintf("Variable `%s` was declared continuous but could not be converted to numeric.", v), call. = FALSE)
}
}
fit_model <- function(selected) {
formula <- make_formula(selected)
if (effect_type == "OR") {
stats::glm(formula, family = stats::binomial("logit"), data = model_data, y = TRUE, x = TRUE, model = TRUE)
} else {
stats::glm(formula, family = stats::poisson("log"), data = model_data, y = TRUE, x = TRUE, model = TRUE)
}
}
covariance_for <- function(fit) {
if (effect_type == "OR") tryCatch(stats::vcov(fit), error = function(e) NULL) else robust_vcov_poisson(fit)
}
selected_terms <- function(fit) intersect(attr(stats::terms(fit), "term.labels"), vars$variable)
# Whole-variable global LRT p-values.
global_p <- function(fit) {
if (!length(attr(stats::terms(fit), "term.labels"))) return(stats::setNames(numeric(), character()))
z <- tryCatch(suppressWarnings(stats::drop1(fit, test = "LRT")), error = function(e) NULL)
if (is.null(z) || !"Pr(>Chi)" %in% colnames(z)) return(stats::setNames(numeric(), character()))
keep <- rownames(z) != "<none>"
ans <- z[keep, "Pr(>Chi)"]
names(ans) <- rownames(z)[keep]
ans
}
# ----------------------------- selectors -----------------------------
select_step <- function(direction) {
k <- if (criterion == "AIC") 2 else log(nrow(model_data))
lower_fit <- fit_model(forced)
upper_formula <- make_formula(vars$variable)
lower_formula <- make_formula(forced)
fit <- if (direction == "forward") {
tryCatch(stats::step(lower_fit, scope = list(lower = lower_formula, upper = upper_formula), direction = "forward", trace = 0, k = k), error = function(e) lower_fit)
} else {
full_fit <- fit_model(vars$variable)
tryCatch(stats::step(full_fit, scope = list(lower = lower_formula, upper = upper_formula), direction = "backward", trace = 0, k = k), error = function(e) full_fit)
}
unique(c(forced, selected_terms(fit)))
}
select_purposeful <- function() {
candidates <- setdiff(vars$variable, forced)
current <- forced
for (v in candidates) {
fit <- tryCatch(fit_model(unique(c(forced, v))), error = function(e) NULL)
if (is.null(fit)) next
p <- global_p(fit)[v]
if (length(p) && is.finite(p) && p < entry) current <- unique(c(current, v))
}
if (!length(current)) return(character())
protected <- forced
repeat {
fit <- tryCatch(fit_model(current), error = function(e) NULL)
if (is.null(fit)) break
gp <- global_p(fit)
removable <- setdiff(intersect(names(gp), current), protected)
if (!length(removable)) break
pv <- gp[removable]
pv[!is.finite(pv)] <- -Inf
worst <- removable[which.max(pv)]
if (!length(worst) || max(pv) <= stay) break
reduced_vars <- setdiff(current, worst)
reduced_fit <- tryCatch(fit_model(reduced_vars), error = function(e) NULL)
if (is.null(reduced_fit)) break
b1 <- stats::coef(fit)
b0 <- stats::coef(reduced_fit)
common <- setdiff(intersect(names(b1), names(b0)), "(Intercept)")
change <- if (!length(common)) 0 else max(abs(b0[common] - b1[common]) / pmax(abs(b1[common]), 1e-8), na.rm = TRUE)
if (is.finite(change) && change >= confounding) protected <- unique(c(protected, worst)) else current <- reduced_vars
if (setequal(current, protected)) break
}
unique(current)
}
select_lasso <- function() {
if (!requireNamespace("glmnet", quietly = TRUE))
stop("Method 'lasso' requires package `glmnet`; install it with install.packages('glmnet').", call. = FALSE)
full_fit <- fit_model(vars$variable)
Xfull <- stats::model.matrix(full_fit)
all_cols <- colnames(Xfull)
keep <- all_cols != "(Intercept)"
X <- Xfull[, keep, drop = FALSE]
assign_full <- attr(Xfull, "assign")[keep]
terms_full <- attr(stats::terms(full_fit), "term.labels")
penalty <- rep(1, ncol(X))
if (length(forced)) {
forced_index <- match(forced, terms_full)
penalty[assign_full %in% forced_index] <- 0
}
family_name <- if (effect_type == "OR") "binomial" else "poisson"
cv <- glmnet::cv.glmnet(X, model_data$.outcome, family = family_name, alpha = 1, standardize = TRUE, penalty.factor = penalty)
lambda_value <- if (lasso_lambda == "lambda.min") cv$lambda.min else cv$lambda.1se
co <- as.matrix(stats::coef(cv$glmnet.fit, s = lambda_value))
nonzero <- setdiff(rownames(co)[as.numeric(co[, 1L]) != 0], "(Intercept)")
if (!length(nonzero)) return(forced)
column_positions <- match(nonzero, colnames(X))
term_positions <- unique(assign_full[column_positions[!is.na(column_positions)]])
chosen <- terms_full[term_positions]
unique(c(forced, intersect(chosen, vars$variable)))
}
enumerate_models <- function() {
candidates <- setdiff(vars$variable, forced)
if (length(candidates) > max_subset_vars)
stop(sprintf("Exhaustive selection has %d candidate variables; increase `max_subset_vars` or use another method.", length(candidates)), call. = FALSE)
total <- 2^length(candidates)
if (!is.finite(total) || total > max_subset_models)
stop(sprintf("Exhaustive selection requires %.0f models, exceeding `max_subset_models = %.0f`.", total, max_subset_models), call. = FALSE)
sets <- vector("list", total)
fits <- vector("list", total)
scores <- rep(Inf, total)
for (m in 0:(total - 1L)) {
bits <- if (!length(candidates)) logical() else as.logical(intToBits(m)[seq_along(candidates)])
selected <- unique(c(forced, candidates[bits]))
fit <- tryCatch(fit_model(selected), error = function(e) NULL)
sets[[m + 1L]] <- selected
fits[[m + 1L]] <- fit
if (!is.null(fit)) scores[m + 1L] <- if (criterion == "AIC") stats::AIC(fit) else stats::BIC(fit)
}
valid <- is.finite(scores)
list(sets = sets[valid], fits = fits[valid], scores = scores[valid])
}
select_best <- function() {
z <- enumerate_models()
if (!length(z$scores)) stop("No valid subset model could be fitted.", call. = FALSE)
z$sets[[which.min(z$scores)]]
}
select_bma <- function() {
z <- enumerate_models()
if (!length(z$scores)) stop("No valid subset model could be fitted.", call. = FALSE)
bic <- vapply(z$fits, stats::BIC, numeric(1))
delta <- bic - min(bic)
w <- exp(-0.5 * delta)
w <- w / sum(w)
pip <- stats::setNames(rep(0, nrow(vars)), vars$variable)
for (i in seq_along(z$sets)) pip[z$sets[[i]]] <- pip[z$sets[[i]]] + w[i]
chosen <- names(pip)[pip >= bma_pip]
chosen <- unique(c(forced, chosen))
attr(chosen, "pip") <- pip
chosen
}
selected_by_method <- list()
fitted_by_method <- list()
selection_details <- list()
for (m in methods) {
chosen <- switch(m,
full = vars$variable,
forward = select_step("forward"),
backward = select_step("backward"),
purposeful = select_purposeful(),
lasso = select_lasso(),
bma = select_bma(),
best = select_best()
)
details <- attributes(chosen)
chosen <- vars$variable[vars$variable %in% unique(as.character(chosen))]
chosen <- vars$variable[vars$variable %in% unique(c(forced, chosen))]
fit <- tryCatch(fit_model(chosen), error = function(e) stop(sprintf("Refitting method '%s' failed: %s", m, e$message), call. = FALSE))
selected_by_method[[m]] <- chosen
fitted_by_method[[m]] <- fit
selection_details[[m]] <- details
}
# ----------------------------- extract results -----------------------------
extract_fit <- function(fit, selected) {
beta <- stats::coef(fit)
covariance <- covariance_for(fit)
if (is.null(covariance)) stop("The covariance matrix could not be estimated.", call. = FALSE)
se <- sqrt(diag(covariance))
z <- beta / se
p <- 2 * stats::pnorm(abs(z), lower.tail = FALSE)
estimate <- exp(beta)
lower <- exp(beta - 1.96 * se)
upper <- exp(beta + 1.96 * se)
coef_table <- data.frame(term = names(beta), beta = unname(beta), se = unname(se), estimate = unname(estimate), lower = unname(lower), upper = unname(upper), p = unname(p), stringsAsFactors = FALSE)
gp <- global_p(fit)
list(coef = coef_table, global_p = gp)
}
extracted <- Map(extract_fit, fitted_by_method, selected_by_method)
# Reliable coefficient mapping through model.matrix assign and contrasts.
effect_for <- function(method, variable, level = NULL) {
fit <- fitted_by_method[[method]]
ext <- extracted[[method]]
if (!variable %in% selected_by_method[[method]]) return(NULL)
term_labels <- attr(stats::terms(fit), "term.labels")
term_index <- match(variable, term_labels)
if (is.na(term_index)) return(NULL)
mm <- stats::model.matrix(fit)
assign <- attr(mm, "assign")
coef_names <- colnames(mm)[assign == term_index]
coef_names <- intersect(coef_names, ext$coef$term)
meta_index <- match(variable, meta$variable)
if (meta$type[meta_index] != "categorical") {
if (!length(coef_names)) return(NULL)
return(ext$coef[match(coef_names[1L], ext$coef$term), , drop = FALSE])
}
levels_v <- levels(model_data[[variable]])
if (is.null(level) || identical(as.character(level), levels_v[1L])) return("REF")
nonref <- levels_v[-1L]
pos <- match(as.character(level), nonref)
if (is.na(pos) || pos > length(coef_names)) return(NULL)
ext$coef[match(coef_names[pos], ext$coef$term), , drop = FALSE]
}
# Model diagnostics.
auc_value <- function(y, score) {
if (length(unique(y)) != 2L) return(NA_real_)
n1 <- sum(y == 1L); n0 <- sum(y == 0L)
if (!n1 || !n0) return(NA_real_)
r <- rank(score, ties.method = "average")
(sum(r[y == 1L]) - n1 * (n1 + 1) / 2) / (n1 * n0)
}
hosmer_p <- function(fit, groups = 10L) {
if (effect_type != "OR") return(NA_real_)
y0 <- fit$y; pr <- stats::fitted(fit)
if (length(y0) < 20L || length(unique(pr)) < 3L) return(NA_real_)
groups <- min(groups, max(2L, floor(length(y0) / 5L)))
br <- unique(stats::quantile(pr, seq(0, 1, length.out = groups + 1L), names = FALSE))
if (length(br) < 3L) return(NA_real_)
g <- cut(pr, breaks = br, include.lowest = TRUE, labels = FALSE)
obs <- as.numeric(rowsum(y0, g)); exp <- as.numeric(rowsum(pr, g)); nn <- as.numeric(table(g))
stat <- sum((obs - exp)^2 / pmax(exp, .Machine$double.eps) + ((nn - obs) - (nn - exp))^2 / pmax(nn - exp, .Machine$double.eps))
df <- length(nn) - 2L
if (df <= 0L) NA_real_ else stats::pchisq(stat, df, lower.tail = FALSE)
}
vif_range <- function(fit) {
X <- stats::model.matrix(fit)
X <- X[, colnames(X) != "(Intercept)", drop = FALSE]
if (ncol(X) <= 1L) return(c(1, 1))
ok <- apply(X, 2L, function(z) is.finite(stats::sd(z)) && stats::sd(z) > 0)
X <- X[, ok, drop = FALSE]
if (ncol(X) <= 1L) return(c(1, 1))
corx <- suppressWarnings(stats::cor(X))
inv <- tryCatch(solve(corx), error = function(e) tryCatch(qr.solve(corx), error = function(e2) NULL))
if (is.null(inv)) return(c(NA_real_, NA_real_))
values <- diag(inv); values <- values[is.finite(values) & values >= 1]
if (!length(values)) c(NA_real_, NA_real_) else range(values)
}
diagnostics_for <- function(fit, selected) {
null_fit <- tryCatch(stats::update(fit, . ~ 1), error = function(e) NULL)
ll <- as.numeric(stats::logLik(fit)); ll0 <- if (is.null(null_fit)) NA_real_ else as.numeric(stats::logLik(null_fit))
n <- stats::nobs(fit)
mcf <- if (is.finite(ll0) && ll0 != 0) 1 - ll / ll0 else NA_real_
nag <- NA_real_
if (is.finite(ll) && is.finite(ll0)) {
num <- 1 - exp((2 / n) * (ll0 - ll)); den <- 1 - exp((2 / n) * ll0)
if (is.finite(den) && den != 0) nag <- num / den
}
pearson <- sum(stats::residuals(fit, type = "pearson")^2, na.rm = TRUE)
pearson_p <- if (stats::df.residual(fit) > 0L) stats::pchisq(pearson, stats::df.residual(fit), lower.tail = FALSE) else NA_real_
vr <- vif_range(fit)
c(N = n, Events = sum(fit$y == 1L), Variables = length(selected), Parameters = length(stats::coef(fit)),
AIC = stats::AIC(fit), BIC = stats::BIC(fit), McFadden = mcf, Nagelkerke = nag,
Hosmer = hosmer_p(fit), Pearson = pearson_p,
AUC = if (effect_type == "OR") auc_value(fit$y, stats::fitted(fit)) else NA_real_,
VIF_min = vr[1L], VIF_max = vr[2L])
}
diagnostics <- Map(diagnostics_for, fitted_by_method, selected_by_method)
# ----------------------------- rows -----------------------------
rows <- list()
for (i in seq_len(nrow(meta))) {
v <- meta$variable[i]
label <- display_label(get_label(data[[v]], v), v)
if (meta$type[i] == "categorical") {
lev <- levels(model_data[[v]])
if (v %in% reverse_rows) lev <- rev(lev)
rows[[length(rows) + 1L]] <- list(variable = v, label = label, item = "", type = "header")
for (lv in lev) rows[[length(rows) + 1L]] <- list(variable = v, label = label, item = as.character(lv), type = "level")
} else {
rows[[length(rows) + 1L]] <- list(variable = v, label = label, item = "", type = "continuous")
}
}
method_labels <- c(full = "Full", forward = "Forward", backward = "Backward", purposeful = "Purposeful", lasso = "LASSO + refit", bma = "BMA + refit", best = "Best subsets + refit")
cell_effect_html <- function(method, row) {
if (row$type == "header") return("")
z <- effect_for(method, row$variable, if (row$type == "level") row$item else NULL)
if (is.null(z)) return("")
if (is.character(z) && identical(z, "REF")) return("Ref")
escape_html(format_effect(z$estimate[1L], z$lower[1L], z$upper[1L]))
}
cell_p_html <- function(method, row) {
# Global p-values are placed on the variable header row, in the same p
# column used for coefficient p-values. No extra global-p column is added.
if (row$type == "header") {
if (!isTRUE(global) || !row$variable %in% selected_by_method[[method]]) return("")
gp <- extracted[[method]]$global_p[row$variable]
if (!length(gp) || !is.finite(gp)) return("")
return(paste0("Global p = ", format_p_html(gp)))
}
if (row$type == "continuous" && isTRUE(global) && !isTRUE(pvalue)) {
gp <- extracted[[method]]$global_p[row$variable]
if (length(gp) && is.finite(gp)) return(format_p_html(gp))
}
if (!isTRUE(pvalue)) return("")
z <- effect_for(method, row$variable, if (row$type == "level") row$item else NULL)
if (is.null(z) || (is.character(z) && identical(z, "REF"))) return("")
format_p_html(z$p[1L])
}
body <- character(length(rows))
for (i in seq_along(rows)) {
r <- rows[[i]]
characteristic <- if (r$type %in% c("header", "continuous")) {
paste0("<span class=\"variable-name\">", r$label, "</span>")
} else {
paste0("<span class=\"level-name\">", escape_html(r$item), "</span>")
}
cells <- vapply(methods, function(m) {
effect_cell <- paste0("<td class=\"effect\">", cell_effect_html(m, r), "</td>")
p_cell <- if (isTRUE(pvalue) || isTRUE(global)) paste0("<td class=\"effect-p\">", cell_p_html(m, r), "</td>") else ""
paste0(effect_cell, p_cell)
}, character(1))
cls <- paste0("row-", r$type, if (r$type %in% c("header", "continuous")) " variable-start" else "")
body[i] <- paste0("<tr class=\"", cls, "\"><td>", characteristic, "</td>", paste(cells, collapse = ""), "</tr>")
}
# Diagnostics as final table rows.
diagnostic_rows <- c("N", "Events", "Variables", "Parameters", "AIC", "BIC", "McFadden R2", "Nagelkerke R2", "Hosmer-Lemeshow p", "Pearson GOF p", "AUC", "VIF range")
diag_key <- c("N", "Events", "Variables", "Parameters", "AIC", "BIC", "McFadden", "Nagelkerke", "Hosmer", "Pearson", "AUC", "VIF")
diag_body <- character(length(diagnostic_rows))
for (i in seq_along(diagnostic_rows)) {
cells <- vapply(methods, function(m) {
d <- diagnostics[[m]]
key <- diag_key[i]
val <- switch(key,
N = formatC(d["N"], format = "f", digits = 0),
Events = formatC(d["Events"], format = "f", digits = 0),
Variables = formatC(d["Variables"], format = "f", digits = 0),
Parameters = formatC(d["Parameters"], format = "f", digits = 0),
AIC = format_number(d["AIC"], 2), BIC = format_number(d["BIC"], 2),
McFadden = format_number(d["McFadden"], 3), Nagelkerke = format_number(d["Nagelkerke"], 3),
Hosmer = format_p_html(d["Hosmer"]), Pearson = format_p_html(d["Pearson"]),
AUC = format_number(d["AUC"], 3),
VIF = if (all(is.finite(d[c("VIF_min", "VIF_max")]))) paste0(format_number(d["VIF_min"], 2), " - ", format_number(d["VIF_max"], 2)) else ""
)
paste0("<td class=\"diagnostic\">", val, "</td>", if (isTRUE(pvalue) || isTRUE(global)) "<td class=\"diagnostic diagnostic-p\"></td>" else "")
}, character(1))
diag_body[i] <- paste0("<tr class=\"diagnostic-row", if (i == 1L) " diagnostic-start" else "", "\"><td>", escape_html(diagnostic_rows[i]), "</td>", paste(cells, collapse = ""), "</tr>")
}
include_p_column <- isTRUE(pvalue) || isTRUE(global)
top_headings <- paste(vapply(methods,function(m)paste0("<th colspan='",if(include_p_column)2 else 1,"' class='model-header'>",escape_html(method_labels[m]),"</th>"),character(1)),collapse="")
sub_headings <- paste(rep(paste0("<th class='effect-header'>",escape_html(paste0(effect_type," (95% CI)")),"</th>",if(include_p_column)"<th class='p-header'>p-value</th>" else ""),length(methods)),collapse="")
header <- paste0(
"<thead><tr><th rowspan=\"2\">Characteristic</th>", top_headings,
"</tr><tr>", sub_headings, "</tr></thead>"
)
selected_note <- paste(vapply(methods, function(m) paste0(method_labels[m], ": ", if (length(selected_by_method[[m]])) paste(selected_by_method[[m]], collapse = ", ") else "intercept only"), character(1)), collapse = "; ")
notes <- paste0(
"<div class=\"effect-note\">Event = ", escape_html(event_level), ". ",
if (effect_type == "OR") "Logistic regression was used." else "Modified Poisson regression with robust variance was used.",
" Selection procedures selected complete variables and every final model was refitted using the ordinary model. Categorical reference categories are marked Ref; continuous effects are per one-unit increase.</div>",
"<div class=\"model-note\">Selected variables - ", escape_html(selected_note), ". All models used the same complete-case sample (n = ", formatC(nrow(model_data), format = "f", digits = 0), ").</div>"
)
if ("lasso" %in% methods) notes <- paste0(notes, "<div class=\"model-note\">LASSO selection used ", escape_html(lasso_lambda), ".</div>")
if ("bma" %in% methods) notes <- paste0(notes, "<div class=\"model-note\">BMA selection used BIC weights and PIP >= ", format_number(bma_pip, 2), ".</div>")
title_html <- if (is.null(title) || !nzchar(as.character(title)[1L])) "" else paste0("<div class=\"table-title\">", escape_html(as.character(title)[1L]), "</div>")
table_block <- paste0("<section class=\"r4vn-table\">", title_html, "<table>", header, "<tbody>", paste(body, collapse = ""), paste(diag_body, collapse = ""), "</tbody></table>", notes, "</section>")
css <- switch(template,
journal = "body{font-family:'Times New Roman',Times,serif;background:#fff;color:#111;margin:18px}.table-title{font-size:18px;font-weight:700;margin:0 0 8px}table{border-collapse:collapse;width:auto;min-width:760px;border-top:2px solid #111;border-bottom:2px solid #111}th{padding:5px 10px;text-align:right;border-bottom:1.5px solid #111;font-weight:700;white-space:nowrap;background:#fff}th:first-child{text-align:center;min-width:260px}td{padding:4px 10px;vertical-align:top;border:0}tr.variable-start td{border-top:1px solid #aaa}",
clean = "body{font-family:Arial,Helvetica,sans-serif;background:#fff;color:#111;margin:18px}.table-title{font-size:18px;font-weight:700;margin:0 0 8px}table{border-collapse:collapse;width:auto;min-width:760px;border-top:2px solid #222;border-bottom:2px solid #222}th{padding:7px 11px;text-align:right;border-bottom:1.5px solid #222;font-weight:700;white-space:nowrap;background:#f3f3f3}th:first-child{text-align:center;min-width:260px}td{padding:5px 11px;vertical-align:top;border-bottom:1px solid #ddd}tr.variable-start td{border-top:1px solid #999}",
minimal = "body{font-family:Arial,Helvetica,sans-serif;background:#fff;color:#111;margin:18px}.table-title{font-size:17px;font-weight:700;margin:0 0 7px}table{border-collapse:collapse;width:auto;min-width:720px;border-top:1.5px solid #222;border-bottom:1.5px solid #222}th{padding:5px 9px;text-align:right;border-bottom:1px solid #555;font-weight:700;white-space:nowrap;background:#fff}th:first-child{text-align:center;min-width:240px}td{padding:4px 9px;vertical-align:top;border:0}tr.variable-start td{border-top:1px solid #ddd}"
)
common_css <- ".table-wrapper{display:inline-block;max-width:100%;overflow-x:auto}.variable-name{font-weight:700}.variable-code{font-family:Consolas,monospace;font-size:.78em;color:#666;font-weight:400}.level-name{display:inline-block;padding-left:22px;white-space:nowrap}.model-header{text-align:center;font-weight:700}.effect-header{text-align:center;font-weight:700}.p-header{text-align:center;width:72px}.effect{text-align:center;white-space:nowrap}.effect-p{text-align:center;white-space:nowrap;width:72px}.diagnostic{text-align:center}.diagnostic-start td{border-top:1.5px solid #555;padding-top:7px}.diagnostic-row td:first-child{font-style:italic}.effect-note,.model-note{font-size:12px;color:#333;margin-top:6px;line-height:1.35}.table-separator{height:24px}"
blocks <- table_block
if (inherits(append, "r4vn_tabmulti")) blocks <- c(append$blocks, table_block)
document <- paste0("<!DOCTYPE html><html><head><meta charset=\"UTF-8\"><meta name=\"viewport\" content=\"width=device-width,initial-scale=1\"><style>", css, common_css, "</style></head><body><div class=\"table-wrapper\">", paste(blocks, collapse = "<div class=\"table-separator\"></div>"), "</div></body></html>")
if (is.null(file)) file <- tempfile(pattern = "r4vn-tabmulti-", fileext = ".html")
if (!is.character(file) || length(file) != 1L || !nzchar(file)) stop("`file` must be a single valid file path.", call. = FALSE)
if (is.character(append) && length(append) == 1L && file.exists(append)) {
old <- paste(readLines(append, warn = FALSE, encoding = "UTF-8"), collapse = "\n")
if (grepl("</body>", old, fixed = TRUE)) {
document <- sub("</body>", paste0("<div class=\"table-separator\"></div>", table_block, "</body>"), old, fixed = TRUE)
file <- append
}
}
writeLines(enc2utf8(document), file, useBytes = TRUE)
# Build a flat data frame for Word/Excel directly from the already
# calculated display cells. This is local by design; no separate tabledata()
# helper or HTML parsing is required.
clean_export_text <- function(x) {
x <- as.character(x)
x <- gsub("<br\\s*/?>", " ", x, ignore.case = TRUE)
x <- gsub("<[^>]+>", "", x)
x <- gsub("<", "<", x, fixed = TRUE)
x <- gsub(">", ">", x, fixed = TRUE)
x <- gsub(""", "\"", x, fixed = TRUE)
x <- gsub("'", "'", x, fixed = TRUE)
x <- gsub(" ", " ", x, fixed = TRUE)
x <- gsub("&", "&", x, fixed = TRUE)
x <- gsub("[\r\n\t]+", " ", x)
x <- gsub("\\s+", " ", x)
trimws(x)
}
export_method_titles <- unname(method_labels[methods])
export_method_titles[is.na(export_method_titles) | !nzchar(export_method_titles)] <-
methods[is.na(export_method_titles) | !nzchar(export_method_titles)]
export_column_names <- "Characteristic"
for (i in seq_along(methods)) {
export_column_names <- c(export_column_names, paste0(export_method_titles[i], " ", effect_type, " (95% CI)"))
if (include_p_column) export_column_names <- c(export_column_names, paste0(export_method_titles[i], " p-value"))
}
export_rows <- lapply(rows, function(current) {
characteristic <- if (current$type %in% c("header", "continuous")) {
clean_export_text(current$label)
} else {
paste0(" ", clean_export_text(current$item))
}
output_row <- characteristic
for (i in seq_along(methods)) {
method <- methods[i]
output_row <- c(output_row, clean_export_text(cell_effect_html(method, current)))
if (include_p_column) output_row <- c(output_row, clean_export_text(cell_p_html(method, current)))
}
output_row
})
export_diag_labels <- c("N", "Events", "Variables", "Parameters", "AIC", "BIC",
"McFadden R2", "Nagelkerke R2", "Hosmer-Lemeshow p",
"Pearson GOF p", "AUC", "VIF range")
export_diag_keys <- c("N", "Events", "Variables", "Parameters", "AIC", "BIC",
"McFadden", "Nagelkerke", "Hosmer", "Pearson", "AUC", "VIF")
for (j in seq_along(export_diag_labels)) {
output_row <- export_diag_labels[j]
key <- export_diag_keys[j]
for (i in seq_along(methods)) {
d <- diagnostics[[methods[i]]]
value <- switch(
key,
N = formatC(d["N"], format = "f", digits = 0),
Events = formatC(d["Events"], format = "f", digits = 0),
Variables = formatC(d["Variables"], format = "f", digits = 0),
Parameters = formatC(d["Parameters"], format = "f", digits = 0),
AIC = format_number(d["AIC"], 2),
BIC = format_number(d["BIC"], 2),
McFadden = format_number(d["McFadden"], 3),
Nagelkerke = format_number(d["Nagelkerke"], 3),
Hosmer = format_p_plain(d["Hosmer"]),
Pearson = format_p_plain(d["Pearson"]),
AUC = format_number(d["AUC"], 3),
VIF = if (all(is.finite(d[c("VIF_min", "VIF_max")]))) {
paste0(format_number(d["VIF_min"], 2), " - ", format_number(d["VIF_max"], 2))
} else ""
)
output_row <- c(output_row, clean_export_text(value))
if (include_p_column) output_row <- c(output_row, "")
}
export_rows[[length(export_rows) + 1L]] <- output_row
}
if (length(export_rows)) {
table_df <- as.data.frame(do.call(rbind, export_rows), stringsAsFactors = FALSE, check.names = FALSE)
} else {
table_df <- as.data.frame(matrix(character(), nrow = 0L, ncol = length(export_column_names)),
stringsAsFactors = FALSE, check.names = FALSE)
}
names(table_df) <- make.unique(export_column_names, sep = "_")
rownames(table_df) <- NULL
output <- list(
rows = rows,
data = table_df,
metadata = vars,
by = by_name,
event = event_level,
effect = effect_type,
methods = methods,
global = global,
rvrow = reverse_rows,
selected = selected_by_method,
models = fitted_by_method,
diagnostics = diagnostics,
selection_details = selection_details,
raw = if (isTRUE(raw)) extracted else NULL,
html = document,
table_html = table_block,
blocks = blocks,
file = normalizePath(file, winslash = "/", mustWork = TRUE),
call = match.call()
)
class(output) <- "r4vn_tabmulti"
if (isTRUE(show)) {
viewer <- getOption("viewer")
if (is.function(viewer)) viewer(output$file) else utils::browseURL(output$file)
}
invisible(output)
}
#' Print or Reopen an R4VN Model-Comparison Table
#'
#' Opens the HTML file stored in an \code{r4vn_tabmulti} object.
#'
#' @param x An object created by \code{tabmulti()}.
#' @param ... Additional arguments currently ignored.
#' @return The input object, invisibly.
#' @keywords internal
#' @method print r4vn_tabmulti
#' @export
print.r4vn_tabmulti <- function(x, ...) {
viewer <- getOption("viewer")
if (is.function(viewer)) viewer(x$file) else utils::browseURL(x$file)
invisible(x)
}
# Version marker used to verify that the current source file has been loaded.
attr(tabmulti, "r4vn_version") <- "tabmulti-1.1.0-2026-08-03"
Try the R4VN package in your browser
Any scripts or data that you put into this service are public.
R4VN documentation built on Sept. 30, 2026, 5:13 p.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Defines functions print.r4vn_tabmulti tabmulti
Documented in print.r4vn_tabmulti tabmulti
# tabmulti_complete.R
# Rebuilt from the user's tabmulti(1).R and aligned with tab(10).R API/rendering.
# LASSO fix: lasso_lambda is validated explicitly and never passed to match.arg().
#' Compare Multivariable Model-Building Strategies
#'
#' Builds and compares variable-selection strategies for a binary outcome.
#' Each strategy selects complete variables or terms, after which the selected
#' model is refitted using ordinary logistic regression or modified Poisson
#' regression so that conventional OR, RR, or PR estimates, 95% confidence
#' intervals, and p-values can be reported.
#'
#' @usage
#' tabmulti(data = NULL, vars = NULL, by = NULL,
#' methods = c("full", "forward", "backward", "purposeful"),
#' digit = 1, p_digit = 3, effect_digit = 2, global = FALSE,
#' pvalue = TRUE, rvrow = NULL, bold_p = TRUE, p_bold = 0.05,
#' or = FALSE, rr = FALSE, pr = FALSE, event = NULL,
#' criterion = c("AIC", "BIC"), force = NULL, entry = 0.20,
#' stay = 0.05, confounding = 0.10, lasso_lambda = "lambda.1se",
#' bma_pip = 0.50, max_subset_vars = 15L, max_subset_models = 100000L,
#' template = c("journal", "clean", "minimal"), append = NULL,
#' file = NULL, raw = FALSE, name = FALSE, title = NULL, show = TRUE)
#'
#' @param data Optional data frame. When omitted or \code{NULL}, the active
#' data frame set by \code{usedf()} or \code{opendata(..., active = TRUE)}
#' is used.
#' @param vars A variable specification created by \code{vars()}.
#' @param by Binary outcome supplied without quotation marks.
#' @param methods Model strategies: \code{"full"}, \code{"forward"},
#' \code{"backward"}, \code{"purposeful"}, \code{"lasso"},
#' \code{"bma"}, or \code{"best"}.
#' @param digit Retained for API consistency with \code{tab()}.
#' @param p_digit Number of decimal places for p-values.
#' @param effect_digit Number of decimal places for effect estimates and
#' confidence limits.
#' @param global Logical. Display global likelihood-ratio p-values.
#' @param pvalue Logical. Display coefficient p-value columns.
#' @param rvrow Categorical variables whose displayed level order should be
#' reversed. This does not change model reference categories.
#' @param bold_p Logical. Bold p-values smaller than \code{p_bold}.
#' @param p_bold Threshold used when \code{bold_p = TRUE}.
#' @param or Logical. Report odds ratios from logistic regression.
#' @param rr Logical. Report risk ratios from modified Poisson regression.
#' @param pr Logical. Report prevalence ratios from modified Poisson regression.
#' Exactly one of \code{or}, \code{rr}, and \code{pr} must be
#' \code{TRUE}.
#' @param event Event level. The last observed outcome level is used when omitted.
#' @param criterion Selection criterion, \code{"AIC"} or \code{"BIC"}.
#' @param force Variables forced into every selected model. Accepts
#' \code{vars(...)}, \code{c(...)}, a character vector, \code{TRUE}, or
#' \code{"ALL"}.
#' @param entry Univariate entry threshold for purposeful selection.
#' @param stay Multivariable retention threshold for purposeful selection.
#' @param confounding Relative coefficient-change threshold for identifying a
#' confounder during purposeful selection.
#' @param lasso_lambda Either \code{"lambda.min"} or \code{"lambda.1se"}.
#' @param bma_pip Posterior inclusion-probability threshold used by the
#' BIC-weighted BMA strategy.
#' @param max_subset_vars Maximum number of candidate variables for exhaustive
#' subset methods.
#' @param max_subset_models Maximum number of subset models to evaluate.
#' @param template HTML style: \code{"journal"}, \code{"clean"}, or
#' \code{"minimal"}.
#' @param append Optional previous \code{r4vn_tabmulti} object or existing HTML path.
#' @param file Optional output HTML path.
#' @param raw Logical. Retain unformatted coefficients and selection details.
#' @param name Logical. Display original variable names beside labels.
#' @param title Optional table title.
#' @param show Logical. Open the HTML table in the Viewer or browser.
#'
#' @details
#' Available strategies are:
#' \itemize{
#' \item \code{full}: include every candidate variable;
#' \item \code{forward}: forward stepwise selection using AIC or BIC;
#' \item \code{backward}: backward stepwise selection from the full model;
#' \item \code{purposeful}: univariate screening followed by significance
#' and confounding assessment;
#' \item \code{lasso}: selection with \code{glmnet::cv.glmnet()}, followed
#' by ordinary-model refitting; requires the suggested package \pkg{glmnet};
#' \item \code{bma}: BIC-weighted subset averaging and inclusion-probability
#' thresholding;
#' \item \code{best}: select the subset with the smallest AIC or BIC.
#' }
#'
#' All strategies use the same complete-case sample. Diagnostic rows include
#' sample size, events, number of variables and parameters, AIC, BIC,
#' pseudo-R-squared measures, goodness-of-fit tests, AUC where applicable, and
#' the coefficient-level VIF range. Exhaustive methods grow exponentially with
#' the number of candidate variables.
#'
#' @return Invisibly returns an object of class \code{r4vn_tabmulti}. Important
#' components include \code{data}, \code{selected}, \code{models},
#' \code{diagnostics}, \code{file}, \code{html}, and
#' \code{table_html}.
#'
#' @seealso \code{\link{vars}}, \code{\link{tab}}, and
#' \code{\link{tabexport}}.
#' @family R4VN tables
#'
#' @examples
#' set.seed(2026)
#' n <- 180
#' dat <- data.frame(
#' age = round(rnorm(n, 45, 12)),
#' sex = factor(sample(c("Female", "Male"), n, TRUE)),
#' bmi = round(rnorm(n, 23, 3), 1),
#' smoking = factor(sample(c("No", "Yes"), n, TRUE,
#' prob = c(0.70, 0.30))),
#' education = factor(sample(c("Primary", "Secondary", "College"),
#' n, TRUE))
#' )
#' lp <- -3.1 + 0.045 * dat$age + 0.11 * (dat$bmi - 23) +
#' 0.45 * (dat$sex == "Male") + 0.70 * (dat$smoking == "Yes")
#' dat$hypertension <- factor(
#' rbinom(n, 1, plogis(lp)),
#' levels = c(0, 1), labels = c("No", "Yes")
#' )
#'
#' models <- tabmulti(
#' dat,
#' vars = vars(c.age, b2.sex, c.bmi, b2.smoking, b2.education),
#' by = hypertension,
#' methods = c("full", "backward"),
#' or = TRUE,
#' event = "Yes",
#' criterion = "AIC",
#' global = TRUE,
#' show = FALSE
#' )
#' models$selected
#' models$diagnostics
#'
#' \donttest{
#' # LASSO requires the suggested package glmnet.
#' if (requireNamespace("glmnet", quietly = TRUE)) {
#' models_lasso <- tabmulti(
#' dat,
#' vars = vars(c.age, b2.sex, c.bmi, b2.smoking, b2.education),
#' by = hypertension,
#' methods = c("full", "lasso"),
#' or = TRUE,
#' event = "Yes",
#' show = FALSE
#' )
#' }
#' }
#'
#' # Extended usage examples
#' \donttest{
#' set.seed(2026)
#' n <- 400
#' d <- data.frame(
#' sex = factor(sample(c("Female", "Male"), n, TRUE)),
#' age = rnorm(n, 45, 12),
#' bmi = rnorm(n, 23, 3),
#' smoking = factor(sample(c("No", "Yes"), n, TRUE)),
#' education = factor(sample(c("Primary", "Secondary", "College"), n, TRUE))
#' )
#' lp <- -3.2 + 0.045 * d$age + 0.10 * (d$bmi - 23) +
#' 0.45 * (d$sex == "Male") + 0.70 * (d$smoking == "Yes")
#' d$outcome <- factor(rbinom(n, 1, plogis(lp)),
#' levels = 0:1, labels = c("No", "Yes"))
#'
#' # Full multivariable logistic model
#' m1 <- tabmulti(
#' d,
#' vars = vars(b2.sex, c.age, c.bmi, b2.smoking, b2.education),
#' by = outcome,
#' methods = "full",
#' or = TRUE,
#' event = "Yes",
#' show = FALSE
#' )
#'
#' # Compare several model-building strategies in one table
#' m2 <- tabmulti(
#' d,
#' vars = vars(b2.sex, c.age, c.bmi, b2.smoking, b2.education),
#' by = outcome,
#' methods = c("full", "forward", "backward", "purposeful"),
#' or = TRUE,
#' event = "Yes",
#' criterion = "AIC",
#' global = TRUE,
#' show = FALSE
#' )
#' m2$selected
#' m2$diagnostics
#'
#' # Force variables into every selected model and tune purposeful selection
#' tabmulti(
#' d,
#' vars = vars(b2.sex, c.age, c.bmi, b2.smoking, b2.education),
#' by = outcome,
#' methods = "purposeful",
#' force = vars(age, sex),
#' entry = 0.20, stay = 0.05, confounding = 0.10,
#' or = TRUE, event = "Yes", show = FALSE
#' )
#'
#' # Modified Poisson models for RR or PR
#' tabmulti(d, vars = vars(b2.sex, c.age, c.bmi, b2.smoking),
#' by = outcome, methods = "full", rr = TRUE,
#' event = "Yes", show = FALSE)
#' tabmulti(d, vars = vars(b2.sex, c.age, c.bmi, b2.smoking),
#' by = outcome, methods = "full", pr = TRUE,
#' event = "Yes", show = FALSE)
#'
#' # Display and output controls
#' tabmulti(d, vars = vars(b2.sex, c.age, c.bmi, b2.smoking),
#' by = outcome, methods = c("full", "backward"),
#' or = TRUE, event = "Yes", rvrow = vars(smoking),
#' pvalue = TRUE, bold_p = TRUE, p_bold = 0.05,
#' template = "clean", raw = TRUE, name = TRUE,
#' title = "Model-building comparison", show = FALSE)
#'
#' # Active-data syntax
#' usedf(d)
#' tabmulti(vars = vars(b2.sex, c.age, c.bmi, b2.smoking),
#' by = outcome, methods = "full", or = TRUE,
#' event = "Yes", show = FALSE)
#'
#' # LASSO requires glmnet; BMA/best are exhaustive and suit fewer candidates
#' if (requireNamespace("glmnet", quietly = TRUE)) {
#' tabmulti(d, vars = vars(b2.sex, c.age, c.bmi, b2.smoking),
#' by = outcome, methods = c("full", "lasso"), or = TRUE,
#' event = "Yes", lasso_lambda = "lambda.1se", show = FALSE)
#' }
#' tabmulti(d, vars = vars(b2.sex, c.age, c.bmi, b2.smoking),
#' by = outcome, methods = c("best", "bma"), or = TRUE,
#' event = "Yes", criterion = "BIC", bma_pip = 0.50,
#' max_subset_vars = 10, show = FALSE)
#' }
#' @export
tabmulti <- function(data = NULL, vars = NULL, by = NULL,
methods = c("full", "forward", "backward", "purposeful"),
digit = 1, p_digit = 3, effect_digit = 2,
global = FALSE, pvalue = TRUE, rvrow = NULL,
bold_p = TRUE, p_bold = 0.05,
or = FALSE, rr = FALSE, pr = FALSE, event = NULL,
criterion = c("AIC", "BIC"), force = NULL,
entry = 0.20, stay = 0.05, confounding = 0.10,
lasso_lambda = "lambda.1se", bma_pip = 0.50,
max_subset_vars = 15L, max_subset_models = 100000L,
template = c("journal", "clean", "minimal"),
append = NULL, file = NULL, raw = FALSE,
name = FALSE, title = NULL, show = TRUE) {
# Use active data when `data` is omitted.
if (is.null(data)) data <- .r4vn_get_active()
# ----------------------------- validation -----------------------------
if (!is.data.frame(data)) stop("`data` must be a data frame.", call. = FALSE)
if (!inherits(vars, "r4vn_vars")) stop("`vars` must be created using `vars()`.", call. = FALSE)
vars <- .r4vn_resolve_vars_input(
vars, data = data, arg = "vars", default_type = "auto", strict = TRUE
)
validate_integer <- function(x, arg) {
if (!is.numeric(x) || length(x) != 1L || is.na(x) || x < 0 || x != floor(x))
stop(sprintf("`%s` must be a single non-negative integer.", arg), call. = FALSE)
}
validate_flag <- function(x, arg) {
if (!is.logical(x) || length(x) != 1L || is.na(x))
stop(sprintf("`%s` must be TRUE or FALSE.", arg), call. = FALSE)
}
validate_probability <- function(x, arg) {
if (!is.numeric(x) || length(x) != 1L || is.na(x) || x < 0 || x > 1)
stop(sprintf("`%s` must be between 0 and 1.", arg), call. = FALSE)
}
validate_integer(digit, "digit")
validate_integer(p_digit, "p_digit")
validate_integer(effect_digit, "effect_digit")
validate_integer(max_subset_vars, "max_subset_vars")
validate_integer(max_subset_models, "max_subset_models")
for (a in c("global", "pvalue", "bold_p", "or", "rr", "pr", "raw", "name", "show"))
validate_flag(get(a), a)
validate_probability(p_bold, "p_bold")
validate_probability(entry, "entry")
validate_probability(stay, "stay")
validate_probability(confounding, "confounding")
validate_probability(bma_pip, "bma_pip")
template <- match.arg(template)
criterion <- match.arg(criterion)
allowed_methods <- c("full", "forward", "backward", "purposeful", "lasso", "bma", "best")
methods <- unique(tolower(as.character(methods)))
invalid <- setdiff(methods, allowed_methods)
if (!length(methods)) stop("`methods` must contain at least one method.", call. = FALSE)
if (length(invalid)) stop(sprintf("Unsupported methods: %s.", paste(invalid, collapse = ", ")), call. = FALSE)
# Do not use match.arg() here: exact validation avoids the previous ambiguous failure.
lasso_lambda <- as.character(lasso_lambda)[1L]
if (is.na(lasso_lambda) || !lasso_lambda %in% c("lambda.min", "lambda.1se"))
stop("`lasso_lambda` must be exactly 'lambda.min' or 'lambda.1se'.", call. = FALSE)
effects <- c(OR = or, RR = rr, PR = pr)
if (sum(effects) != 1L) stop("Exactly one of `or`, `rr`, or `pr` must be TRUE.", call. = FALSE)
effect_type <- names(effects)[effects][1L]
by_expr <- substitute(by)
if (!is.symbol(by_expr)) stop("`by` must be a single variable name.", call. = FALSE)
by_name <- as.character(by_expr)
duplicated_variables <- unique(vars$variable[duplicated(vars$variable)])
if (length(duplicated_variables))
stop(sprintf("Duplicated variable specification: %s.", paste(duplicated_variables, collapse = ", ")), call. = FALSE)
absent <- setdiff(c(vars$variable, by_name), names(data))
if (length(absent)) stop(sprintf("Variables not found in `data`: %s.", paste(absent, collapse = ", ")), call. = FALSE)
if (by_name %in% vars$variable) stop("The outcome variable cannot also appear in `vars`.", call. = FALSE)
# ----------------------------- helpers -----------------------------
escape_html <- 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)
}
get_label <- function(x, variable) {
z <- attr(x, "label", exact = TRUE)
if (is.null(z) || !length(z) || is.na(z[1L]) || !nzchar(as.character(z[1L]))) variable else as.character(z[1L])
}
display_label <- function(label, variable) {
if (!isTRUE(name)) return(escape_html(label))
paste0(escape_html(label), " <span class=\"variable-code\">[", escape_html(variable), "]</span>")
}
get_levels <- function(x) {
observed <- x[!is.na(x)]
values <- if (is.factor(x)) levels(x) else if (is.logical(x)) c(FALSE, TRUE) else {
z <- unique(observed)
if (is.numeric(z)) sort(z) else z
}
values[values %in% observed]
}
format_number <- function(x, digits = effect_digit) {
if (!length(x) || is.na(x) || !is.finite(x)) return("")
formatC(x, format = "f", digits = digits, big.mark = ",")
}
format_p_plain <- function(x) {
if (!length(x) || is.na(x) || !is.finite(x)) return("")
limit <- 10^(-p_digit)
if (x < limit) paste0("<", formatC(limit, format = "f", digits = p_digit))
else formatC(x, format = "f", digits = p_digit)
}
format_p_html <- function(x) {
z <- escape_html(format_p_plain(x))
if (nzchar(z) && isTRUE(bold_p) && is.finite(x) && x < p_bold) paste0("<strong>", z, "</strong>") else z
}
format_effect <- function(est, low, high) {
if (!all(is.finite(c(est, low, high)))) return("")
paste0(format_number(est), " (", format_number(low), " - ", format_number(high), ")")
}
make_formula <- function(variables) stats::reformulate(variables, response = ".outcome")
robust_vcov_poisson <- function(fit) {
X <- stats::model.matrix(fit)
mu <- stats::fitted(fit)
score <- fit$y - mu
bread <- tryCatch(solve(crossprod(X, X * as.vector(mu))), error = function(e) NULL)
if (is.null(bread)) return(NULL)
meat <- crossprod(X, X * as.vector(score^2))
out <- bread %*% meat %*% bread
dimnames(out) <- list(colnames(X), colnames(X))
out
}
# Parse force with the same accepted styles as tab(adjusted/multi).
parse_force <- function(expression) {
all_names <- vars$variable
if (identical(expression, quote(NULL))) return(character())
if (is.logical(expression) && length(expression) == 1L)
return(if (isTRUE(expression)) all_names else character())
if (is.character(expression)) {
if (length(expression) == 1L && toupper(expression) == "ALL") return(all_names)
return(expression)
}
if (is.symbol(expression)) {
object_name <- as.character(expression)
if (toupper(object_name) == "ALL") return(all_names)
value <- tryCatch(get(object_name, envir = parent.frame(2L)), error = function(e) NULL)
if (inherits(value, "r4vn_vars")) return(value$variable)
if (is.character(value)) return(if (length(value) == 1L && toupper(value) == "ALL") all_names else value)
if (is.logical(value) && length(value) == 1L && isTRUE(value)) return(all_names)
return(object_name)
}
if (is.call(expression) && as.character(expression[[1L]]) == "vars") {
value <- eval(expression, envir = parent.frame(2L))
if (!inherits(value, "r4vn_vars")) stop("`force = vars(...)` is invalid.", call. = FALSE)
return(value$variable)
}
if (is.call(expression) && as.character(expression[[1L]]) == "c") {
items <- as.list(expression)[-1L]
if (all(vapply(items, is.symbol, logical(1)))) return(vapply(items, as.character, character(1)))
value <- eval(expression, envir = parent.frame(2L))
if (is.character(value)) return(value)
}
value <- tryCatch(eval(expression, envir = parent.frame(2L)), error = function(e) NULL)
if (inherits(value, "r4vn_vars")) return(value$variable)
if (is.character(value)) return(value)
stop("`force` must be NULL, TRUE, 'ALL', vars(...), c(...), or a character vector.", call. = FALSE)
}
forced <- unique(parse_force(substitute(force)))
invalid_forced <- setdiff(forced, vars$variable)
if (length(invalid_forced)) stop(sprintf("Forced variables not listed in `vars`: %s.", paste(invalid_forced, collapse = ", ")), call. = FALSE)
# Parse rvrow using the same conventions as tab(). rvrow changes display
# order only; it never changes the reference category used in a model.
parse_rvrow <- function(expression) {
categorical <- vars$variable[vars$type == "categorical"]
if (identical(expression, quote(NULL))) return(character())
if (is.logical(expression) && length(expression) == 1L)
return(if (isTRUE(expression)) categorical else character())
if (is.character(expression)) return(expression)
if (is.symbol(expression)) {
object_name <- as.character(expression)
value <- tryCatch(get(object_name, envir = parent.frame(2L)), error = function(e) NULL)
if (inherits(value, "r4vn_vars")) return(value$variable)
if (is.character(value)) return(value)
if (is.logical(value) && length(value) == 1L) return(if (isTRUE(value)) categorical else character())
return(object_name)
}
if (is.call(expression) && as.character(expression[[1L]]) %in% c("vars", "c")) {
if (as.character(expression[[1L]]) == "vars") {
value <- eval(expression, envir = parent.frame(2L))
if (inherits(value, "r4vn_vars")) return(value$variable)
} else {
items <- as.list(expression)[-1L]
if (all(vapply(items, is.symbol, logical(1)))) return(vapply(items, as.character, character(1)))
value <- tryCatch(eval(expression, envir = parent.frame(2L)), error = function(e) NULL)
if (is.character(value)) return(value)
}
}
value <- tryCatch(eval(expression, envir = parent.frame(2L)), error = function(e) NULL)
if (inherits(value, "r4vn_vars")) return(value$variable)
if (is.character(value)) return(value)
stop("`rvrow` must be NULL, TRUE, vars(...), c(...), a variable name, or a character vector.", call. = FALSE)
}
reverse_rows <- unique(parse_rvrow(substitute(rvrow)))
invalid_reverse <- setdiff(reverse_rows, vars$variable[vars$type == "categorical"])
if (length(invalid_reverse)) stop(sprintf("Invalid categorical variables in `rvrow`: %s.", paste(invalid_reverse, collapse = ", ")), call. = FALSE)
# ----------------------------- outcome/data -----------------------------
y <- data[[by_name]]
original_y_levels <- get_levels(y)
if (length(original_y_levels) != 2L) stop("`by` must have exactly two observed levels.", call. = FALSE)
event_level <- if (is.null(event)) utils::tail(original_y_levels, 1L) else as.character(event)[1L]
if (!event_level %in% as.character(original_y_levels)) stop("`event` is not an observed outcome level.", call. = FALSE)
model_data <- data[c(by_name, vars$variable)]
model_data <- model_data[stats::complete.cases(model_data), , drop = FALSE]
if (!nrow(model_data)) stop("No complete observations remain for the requested variables.", call. = FALSE)
model_data$.outcome <- as.integer(as.character(model_data[[by_name]]) == event_level)
if (length(unique(model_data$.outcome)) != 2L) stop("The complete-case data do not contain both outcome groups.", call. = FALSE)
# Preserve metadata and references exactly as declared by vars().
meta <- vars
for (i in seq_len(nrow(meta))) {
v <- meta$variable[i]
if (meta$type[i] == "categorical") {
lev <- get_levels(data[[v]])
ref_index <- meta$reference_index[i]
if (is.na(ref_index) || ref_index < 1L || ref_index > length(lev))
stop(sprintf("Invalid reference index for `%s`.", v), call. = FALSE)
ref <- as.character(lev[ref_index])
model_data[[v]] <- factor(model_data[[v]], levels = lev)
model_data[[v]] <- stats::relevel(model_data[[v]], ref = ref)
} else {
model_data[[v]] <- suppressWarnings(as.numeric(model_data[[v]]))
if (anyNA(model_data[[v]])) stop(sprintf("Variable `%s` was declared continuous but could not be converted to numeric.", v), call. = FALSE)
}
}
fit_model <- function(selected) {
formula <- make_formula(selected)
if (effect_type == "OR") {
stats::glm(formula, family = stats::binomial("logit"), data = model_data, y = TRUE, x = TRUE, model = TRUE)
} else {
stats::glm(formula, family = stats::poisson("log"), data = model_data, y = TRUE, x = TRUE, model = TRUE)
}
}
covariance_for <- function(fit) {
if (effect_type == "OR") tryCatch(stats::vcov(fit), error = function(e) NULL) else robust_vcov_poisson(fit)
}
selected_terms <- function(fit) intersect(attr(stats::terms(fit), "term.labels"), vars$variable)
# Whole-variable global LRT p-values.
global_p <- function(fit) {
if (!length(attr(stats::terms(fit), "term.labels"))) return(stats::setNames(numeric(), character()))
z <- tryCatch(suppressWarnings(stats::drop1(fit, test = "LRT")), error = function(e) NULL)
if (is.null(z) || !"Pr(>Chi)" %in% colnames(z)) return(stats::setNames(numeric(), character()))
keep <- rownames(z) != "<none>"
ans <- z[keep, "Pr(>Chi)"]
names(ans) <- rownames(z)[keep]
ans
}
# ----------------------------- selectors -----------------------------
select_step <- function(direction) {
k <- if (criterion == "AIC") 2 else log(nrow(model_data))
lower_fit <- fit_model(forced)
upper_formula <- make_formula(vars$variable)
lower_formula <- make_formula(forced)
fit <- if (direction == "forward") {
tryCatch(stats::step(lower_fit, scope = list(lower = lower_formula, upper = upper_formula), direction = "forward", trace = 0, k = k), error = function(e) lower_fit)
} else {
full_fit <- fit_model(vars$variable)
tryCatch(stats::step(full_fit, scope = list(lower = lower_formula, upper = upper_formula), direction = "backward", trace = 0, k = k), error = function(e) full_fit)
}
unique(c(forced, selected_terms(fit)))
}
select_purposeful <- function() {
candidates <- setdiff(vars$variable, forced)
current <- forced
for (v in candidates) {
fit <- tryCatch(fit_model(unique(c(forced, v))), error = function(e) NULL)
if (is.null(fit)) next
p <- global_p(fit)[v]
if (length(p) && is.finite(p) && p < entry) current <- unique(c(current, v))
}
if (!length(current)) return(character())
protected <- forced
repeat {
fit <- tryCatch(fit_model(current), error = function(e) NULL)
if (is.null(fit)) break
gp <- global_p(fit)
removable <- setdiff(intersect(names(gp), current), protected)
if (!length(removable)) break
pv <- gp[removable]
pv[!is.finite(pv)] <- -Inf
worst <- removable[which.max(pv)]
if (!length(worst) || max(pv) <= stay) break
reduced_vars <- setdiff(current, worst)
reduced_fit <- tryCatch(fit_model(reduced_vars), error = function(e) NULL)
if (is.null(reduced_fit)) break
b1 <- stats::coef(fit)
b0 <- stats::coef(reduced_fit)
common <- setdiff(intersect(names(b1), names(b0)), "(Intercept)")
change <- if (!length(common)) 0 else max(abs(b0[common] - b1[common]) / pmax(abs(b1[common]), 1e-8), na.rm = TRUE)
if (is.finite(change) && change >= confounding) protected <- unique(c(protected, worst)) else current <- reduced_vars
if (setequal(current, protected)) break
}
unique(current)
}
select_lasso <- function() {
if (!requireNamespace("glmnet", quietly = TRUE))
stop("Method 'lasso' requires package `glmnet`; install it with install.packages('glmnet').", call. = FALSE)
full_fit <- fit_model(vars$variable)
Xfull <- stats::model.matrix(full_fit)
all_cols <- colnames(Xfull)
keep <- all_cols != "(Intercept)"
X <- Xfull[, keep, drop = FALSE]
assign_full <- attr(Xfull, "assign")[keep]
terms_full <- attr(stats::terms(full_fit), "term.labels")
penalty <- rep(1, ncol(X))
if (length(forced)) {
forced_index <- match(forced, terms_full)
penalty[assign_full %in% forced_index] <- 0
}
family_name <- if (effect_type == "OR") "binomial" else "poisson"
cv <- glmnet::cv.glmnet(X, model_data$.outcome, family = family_name, alpha = 1, standardize = TRUE, penalty.factor = penalty)
lambda_value <- if (lasso_lambda == "lambda.min") cv$lambda.min else cv$lambda.1se
co <- as.matrix(stats::coef(cv$glmnet.fit, s = lambda_value))
nonzero <- setdiff(rownames(co)[as.numeric(co[, 1L]) != 0], "(Intercept)")
if (!length(nonzero)) return(forced)
column_positions <- match(nonzero, colnames(X))
term_positions <- unique(assign_full[column_positions[!is.na(column_positions)]])
chosen <- terms_full[term_positions]
unique(c(forced, intersect(chosen, vars$variable)))
}
enumerate_models <- function() {
candidates <- setdiff(vars$variable, forced)
if (length(candidates) > max_subset_vars)
stop(sprintf("Exhaustive selection has %d candidate variables; increase `max_subset_vars` or use another method.", length(candidates)), call. = FALSE)
total <- 2^length(candidates)
if (!is.finite(total) || total > max_subset_models)
stop(sprintf("Exhaustive selection requires %.0f models, exceeding `max_subset_models = %.0f`.", total, max_subset_models), call. = FALSE)
sets <- vector("list", total)
fits <- vector("list", total)
scores <- rep(Inf, total)
for (m in 0:(total - 1L)) {
bits <- if (!length(candidates)) logical() else as.logical(intToBits(m)[seq_along(candidates)])
selected <- unique(c(forced, candidates[bits]))
fit <- tryCatch(fit_model(selected), error = function(e) NULL)
sets[[m + 1L]] <- selected
fits[[m + 1L]] <- fit
if (!is.null(fit)) scores[m + 1L] <- if (criterion == "AIC") stats::AIC(fit) else stats::BIC(fit)
}
valid <- is.finite(scores)
list(sets = sets[valid], fits = fits[valid], scores = scores[valid])
}
select_best <- function() {
z <- enumerate_models()
if (!length(z$scores)) stop("No valid subset model could be fitted.", call. = FALSE)
z$sets[[which.min(z$scores)]]
}
select_bma <- function() {
z <- enumerate_models()
if (!length(z$scores)) stop("No valid subset model could be fitted.", call. = FALSE)
bic <- vapply(z$fits, stats::BIC, numeric(1))
delta <- bic - min(bic)
w <- exp(-0.5 * delta)
w <- w / sum(w)
pip <- stats::setNames(rep(0, nrow(vars)), vars$variable)
for (i in seq_along(z$sets)) pip[z$sets[[i]]] <- pip[z$sets[[i]]] + w[i]
chosen <- names(pip)[pip >= bma_pip]
chosen <- unique(c(forced, chosen))
attr(chosen, "pip") <- pip
chosen
}
selected_by_method <- list()
fitted_by_method <- list()
selection_details <- list()
for (m in methods) {
chosen <- switch(m,
full = vars$variable,
forward = select_step("forward"),
backward = select_step("backward"),
purposeful = select_purposeful(),
lasso = select_lasso(),
bma = select_bma(),
best = select_best()
)
details <- attributes(chosen)
chosen <- vars$variable[vars$variable %in% unique(as.character(chosen))]
chosen <- vars$variable[vars$variable %in% unique(c(forced, chosen))]
fit <- tryCatch(fit_model(chosen), error = function(e) stop(sprintf("Refitting method '%s' failed: %s", m, e$message), call. = FALSE))
selected_by_method[[m]] <- chosen
fitted_by_method[[m]] <- fit
selection_details[[m]] <- details
}
# ----------------------------- extract results -----------------------------
extract_fit <- function(fit, selected) {
beta <- stats::coef(fit)
covariance <- covariance_for(fit)
if (is.null(covariance)) stop("The covariance matrix could not be estimated.", call. = FALSE)
se <- sqrt(diag(covariance))
z <- beta / se
p <- 2 * stats::pnorm(abs(z), lower.tail = FALSE)
estimate <- exp(beta)
lower <- exp(beta - 1.96 * se)
upper <- exp(beta + 1.96 * se)
coef_table <- data.frame(term = names(beta), beta = unname(beta), se = unname(se), estimate = unname(estimate), lower = unname(lower), upper = unname(upper), p = unname(p), stringsAsFactors = FALSE)
gp <- global_p(fit)
list(coef = coef_table, global_p = gp)
}
extracted <- Map(extract_fit, fitted_by_method, selected_by_method)
# Reliable coefficient mapping through model.matrix assign and contrasts.
effect_for <- function(method, variable, level = NULL) {
fit <- fitted_by_method[[method]]
ext <- extracted[[method]]
if (!variable %in% selected_by_method[[method]]) return(NULL)
term_labels <- attr(stats::terms(fit), "term.labels")
term_index <- match(variable, term_labels)
if (is.na(term_index)) return(NULL)
mm <- stats::model.matrix(fit)
assign <- attr(mm, "assign")
coef_names <- colnames(mm)[assign == term_index]
coef_names <- intersect(coef_names, ext$coef$term)
meta_index <- match(variable, meta$variable)
if (meta$type[meta_index] != "categorical") {
if (!length(coef_names)) return(NULL)
return(ext$coef[match(coef_names[1L], ext$coef$term), , drop = FALSE])
}
levels_v <- levels(model_data[[variable]])
if (is.null(level) || identical(as.character(level), levels_v[1L])) return("REF")
nonref <- levels_v[-1L]
pos <- match(as.character(level), nonref)
if (is.na(pos) || pos > length(coef_names)) return(NULL)
ext$coef[match(coef_names[pos], ext$coef$term), , drop = FALSE]
}
# Model diagnostics.
auc_value <- function(y, score) {
if (length(unique(y)) != 2L) return(NA_real_)
n1 <- sum(y == 1L); n0 <- sum(y == 0L)
if (!n1 || !n0) return(NA_real_)
r <- rank(score, ties.method = "average")
(sum(r[y == 1L]) - n1 * (n1 + 1) / 2) / (n1 * n0)
}
hosmer_p <- function(fit, groups = 10L) {
if (effect_type != "OR") return(NA_real_)
y0 <- fit$y; pr <- stats::fitted(fit)
if (length(y0) < 20L || length(unique(pr)) < 3L) return(NA_real_)
groups <- min(groups, max(2L, floor(length(y0) / 5L)))
br <- unique(stats::quantile(pr, seq(0, 1, length.out = groups + 1L), names = FALSE))
if (length(br) < 3L) return(NA_real_)
g <- cut(pr, breaks = br, include.lowest = TRUE, labels = FALSE)
obs <- as.numeric(rowsum(y0, g)); exp <- as.numeric(rowsum(pr, g)); nn <- as.numeric(table(g))
stat <- sum((obs - exp)^2 / pmax(exp, .Machine$double.eps) + ((nn - obs) - (nn - exp))^2 / pmax(nn - exp, .Machine$double.eps))
df <- length(nn) - 2L
if (df <= 0L) NA_real_ else stats::pchisq(stat, df, lower.tail = FALSE)
}
vif_range <- function(fit) {
X <- stats::model.matrix(fit)
X <- X[, colnames(X) != "(Intercept)", drop = FALSE]
if (ncol(X) <= 1L) return(c(1, 1))
ok <- apply(X, 2L, function(z) is.finite(stats::sd(z)) && stats::sd(z) > 0)
X <- X[, ok, drop = FALSE]
if (ncol(X) <= 1L) return(c(1, 1))
corx <- suppressWarnings(stats::cor(X))
inv <- tryCatch(solve(corx), error = function(e) tryCatch(qr.solve(corx), error = function(e2) NULL))
if (is.null(inv)) return(c(NA_real_, NA_real_))
values <- diag(inv); values <- values[is.finite(values) & values >= 1]
if (!length(values)) c(NA_real_, NA_real_) else range(values)
}
diagnostics_for <- function(fit, selected) {
null_fit <- tryCatch(stats::update(fit, . ~ 1), error = function(e) NULL)
ll <- as.numeric(stats::logLik(fit)); ll0 <- if (is.null(null_fit)) NA_real_ else as.numeric(stats::logLik(null_fit))
n <- stats::nobs(fit)
mcf <- if (is.finite(ll0) && ll0 != 0) 1 - ll / ll0 else NA_real_
nag <- NA_real_
if (is.finite(ll) && is.finite(ll0)) {
num <- 1 - exp((2 / n) * (ll0 - ll)); den <- 1 - exp((2 / n) * ll0)
if (is.finite(den) && den != 0) nag <- num / den
}
pearson <- sum(stats::residuals(fit, type = "pearson")^2, na.rm = TRUE)
pearson_p <- if (stats::df.residual(fit) > 0L) stats::pchisq(pearson, stats::df.residual(fit), lower.tail = FALSE) else NA_real_
vr <- vif_range(fit)
c(N = n, Events = sum(fit$y == 1L), Variables = length(selected), Parameters = length(stats::coef(fit)),
AIC = stats::AIC(fit), BIC = stats::BIC(fit), McFadden = mcf, Nagelkerke = nag,
Hosmer = hosmer_p(fit), Pearson = pearson_p,
AUC = if (effect_type == "OR") auc_value(fit$y, stats::fitted(fit)) else NA_real_,
VIF_min = vr[1L], VIF_max = vr[2L])
}
diagnostics <- Map(diagnostics_for, fitted_by_method, selected_by_method)
# ----------------------------- rows -----------------------------
rows <- list()
for (i in seq_len(nrow(meta))) {
v <- meta$variable[i]
label <- display_label(get_label(data[[v]], v), v)
if (meta$type[i] == "categorical") {
lev <- levels(model_data[[v]])
if (v %in% reverse_rows) lev <- rev(lev)
rows[[length(rows) + 1L]] <- list(variable = v, label = label, item = "", type = "header")
for (lv in lev) rows[[length(rows) + 1L]] <- list(variable = v, label = label, item = as.character(lv), type = "level")
} else {
rows[[length(rows) + 1L]] <- list(variable = v, label = label, item = "", type = "continuous")
}
}
method_labels <- c(full = "Full", forward = "Forward", backward = "Backward", purposeful = "Purposeful", lasso = "LASSO + refit", bma = "BMA + refit", best = "Best subsets + refit")
cell_effect_html <- function(method, row) {
if (row$type == "header") return("")
z <- effect_for(method, row$variable, if (row$type == "level") row$item else NULL)
if (is.null(z)) return("")
if (is.character(z) && identical(z, "REF")) return("Ref")
escape_html(format_effect(z$estimate[1L], z$lower[1L], z$upper[1L]))
}
cell_p_html <- function(method, row) {
# Global p-values are placed on the variable header row, in the same p
# column used for coefficient p-values. No extra global-p column is added.
if (row$type == "header") {
if (!isTRUE(global) || !row$variable %in% selected_by_method[[method]]) return("")
gp <- extracted[[method]]$global_p[row$variable]
if (!length(gp) || !is.finite(gp)) return("")
return(paste0("Global p = ", format_p_html(gp)))
}
if (row$type == "continuous" && isTRUE(global) && !isTRUE(pvalue)) {
gp <- extracted[[method]]$global_p[row$variable]
if (length(gp) && is.finite(gp)) return(format_p_html(gp))
}
if (!isTRUE(pvalue)) return("")
z <- effect_for(method, row$variable, if (row$type == "level") row$item else NULL)
if (is.null(z) || (is.character(z) && identical(z, "REF"))) return("")
format_p_html(z$p[1L])
}
body <- character(length(rows))
for (i in seq_along(rows)) {
r <- rows[[i]]
characteristic <- if (r$type %in% c("header", "continuous")) {
paste0("<span class=\"variable-name\">", r$label, "</span>")
} else {
paste0("<span class=\"level-name\">", escape_html(r$item), "</span>")
}
cells <- vapply(methods, function(m) {
effect_cell <- paste0("<td class=\"effect\">", cell_effect_html(m, r), "</td>")
p_cell <- if (isTRUE(pvalue) || isTRUE(global)) paste0("<td class=\"effect-p\">", cell_p_html(m, r), "</td>") else ""
paste0(effect_cell, p_cell)
}, character(1))
cls <- paste0("row-", r$type, if (r$type %in% c("header", "continuous")) " variable-start" else "")
body[i] <- paste0("<tr class=\"", cls, "\"><td>", characteristic, "</td>", paste(cells, collapse = ""), "</tr>")
}
# Diagnostics as final table rows.
diagnostic_rows <- c("N", "Events", "Variables", "Parameters", "AIC", "BIC", "McFadden R2", "Nagelkerke R2", "Hosmer-Lemeshow p", "Pearson GOF p", "AUC", "VIF range")
diag_key <- c("N", "Events", "Variables", "Parameters", "AIC", "BIC", "McFadden", "Nagelkerke", "Hosmer", "Pearson", "AUC", "VIF")
diag_body <- character(length(diagnostic_rows))
for (i in seq_along(diagnostic_rows)) {
cells <- vapply(methods, function(m) {
d <- diagnostics[[m]]
key <- diag_key[i]
val <- switch(key,
N = formatC(d["N"], format = "f", digits = 0),
Events = formatC(d["Events"], format = "f", digits = 0),
Variables = formatC(d["Variables"], format = "f", digits = 0),
Parameters = formatC(d["Parameters"], format = "f", digits = 0),
AIC = format_number(d["AIC"], 2), BIC = format_number(d["BIC"], 2),
McFadden = format_number(d["McFadden"], 3), Nagelkerke = format_number(d["Nagelkerke"], 3),
Hosmer = format_p_html(d["Hosmer"]), Pearson = format_p_html(d["Pearson"]),
AUC = format_number(d["AUC"], 3),
VIF = if (all(is.finite(d[c("VIF_min", "VIF_max")]))) paste0(format_number(d["VIF_min"], 2), " - ", format_number(d["VIF_max"], 2)) else ""
)
paste0("<td class=\"diagnostic\">", val, "</td>", if (isTRUE(pvalue) || isTRUE(global)) "<td class=\"diagnostic diagnostic-p\"></td>" else "")
}, character(1))
diag_body[i] <- paste0("<tr class=\"diagnostic-row", if (i == 1L) " diagnostic-start" else "", "\"><td>", escape_html(diagnostic_rows[i]), "</td>", paste(cells, collapse = ""), "</tr>")
}
include_p_column <- isTRUE(pvalue) || isTRUE(global)
top_headings <- paste(vapply(methods,function(m)paste0("<th colspan='",if(include_p_column)2 else 1,"' class='model-header'>",escape_html(method_labels[m]),"</th>"),character(1)),collapse="")
sub_headings <- paste(rep(paste0("<th class='effect-header'>",escape_html(paste0(effect_type," (95% CI)")),"</th>",if(include_p_column)"<th class='p-header'>p-value</th>" else ""),length(methods)),collapse="")
header <- paste0(
"<thead><tr><th rowspan=\"2\">Characteristic</th>", top_headings,
"</tr><tr>", sub_headings, "</tr></thead>"
)
selected_note <- paste(vapply(methods, function(m) paste0(method_labels[m], ": ", if (length(selected_by_method[[m]])) paste(selected_by_method[[m]], collapse = ", ") else "intercept only"), character(1)), collapse = "; ")
notes <- paste0(
"<div class=\"effect-note\">Event = ", escape_html(event_level), ". ",
if (effect_type == "OR") "Logistic regression was used." else "Modified Poisson regression with robust variance was used.",
" Selection procedures selected complete variables and every final model was refitted using the ordinary model. Categorical reference categories are marked Ref; continuous effects are per one-unit increase.</div>",
"<div class=\"model-note\">Selected variables - ", escape_html(selected_note), ". All models used the same complete-case sample (n = ", formatC(nrow(model_data), format = "f", digits = 0), ").</div>"
)
if ("lasso" %in% methods) notes <- paste0(notes, "<div class=\"model-note\">LASSO selection used ", escape_html(lasso_lambda), ".</div>")
if ("bma" %in% methods) notes <- paste0(notes, "<div class=\"model-note\">BMA selection used BIC weights and PIP >= ", format_number(bma_pip, 2), ".</div>")
title_html <- if (is.null(title) || !nzchar(as.character(title)[1L])) "" else paste0("<div class=\"table-title\">", escape_html(as.character(title)[1L]), "</div>")
table_block <- paste0("<section class=\"r4vn-table\">", title_html, "<table>", header, "<tbody>", paste(body, collapse = ""), paste(diag_body, collapse = ""), "</tbody></table>", notes, "</section>")
css <- switch(template,
journal = "body{font-family:'Times New Roman',Times,serif;background:#fff;color:#111;margin:18px}.table-title{font-size:18px;font-weight:700;margin:0 0 8px}table{border-collapse:collapse;width:auto;min-width:760px;border-top:2px solid #111;border-bottom:2px solid #111}th{padding:5px 10px;text-align:right;border-bottom:1.5px solid #111;font-weight:700;white-space:nowrap;background:#fff}th:first-child{text-align:center;min-width:260px}td{padding:4px 10px;vertical-align:top;border:0}tr.variable-start td{border-top:1px solid #aaa}",
clean = "body{font-family:Arial,Helvetica,sans-serif;background:#fff;color:#111;margin:18px}.table-title{font-size:18px;font-weight:700;margin:0 0 8px}table{border-collapse:collapse;width:auto;min-width:760px;border-top:2px solid #222;border-bottom:2px solid #222}th{padding:7px 11px;text-align:right;border-bottom:1.5px solid #222;font-weight:700;white-space:nowrap;background:#f3f3f3}th:first-child{text-align:center;min-width:260px}td{padding:5px 11px;vertical-align:top;border-bottom:1px solid #ddd}tr.variable-start td{border-top:1px solid #999}",
minimal = "body{font-family:Arial,Helvetica,sans-serif;background:#fff;color:#111;margin:18px}.table-title{font-size:17px;font-weight:700;margin:0 0 7px}table{border-collapse:collapse;width:auto;min-width:720px;border-top:1.5px solid #222;border-bottom:1.5px solid #222}th{padding:5px 9px;text-align:right;border-bottom:1px solid #555;font-weight:700;white-space:nowrap;background:#fff}th:first-child{text-align:center;min-width:240px}td{padding:4px 9px;vertical-align:top;border:0}tr.variable-start td{border-top:1px solid #ddd}"
)
common_css <- ".table-wrapper{display:inline-block;max-width:100%;overflow-x:auto}.variable-name{font-weight:700}.variable-code{font-family:Consolas,monospace;font-size:.78em;color:#666;font-weight:400}.level-name{display:inline-block;padding-left:22px;white-space:nowrap}.model-header{text-align:center;font-weight:700}.effect-header{text-align:center;font-weight:700}.p-header{text-align:center;width:72px}.effect{text-align:center;white-space:nowrap}.effect-p{text-align:center;white-space:nowrap;width:72px}.diagnostic{text-align:center}.diagnostic-start td{border-top:1.5px solid #555;padding-top:7px}.diagnostic-row td:first-child{font-style:italic}.effect-note,.model-note{font-size:12px;color:#333;margin-top:6px;line-height:1.35}.table-separator{height:24px}"
blocks <- table_block
if (inherits(append, "r4vn_tabmulti")) blocks <- c(append$blocks, table_block)
document <- paste0("<!DOCTYPE html><html><head><meta charset=\"UTF-8\"><meta name=\"viewport\" content=\"width=device-width,initial-scale=1\"><style>", css, common_css, "</style></head><body><div class=\"table-wrapper\">", paste(blocks, collapse = "<div class=\"table-separator\"></div>"), "</div></body></html>")
if (is.null(file)) file <- tempfile(pattern = "r4vn-tabmulti-", fileext = ".html")
if (!is.character(file) || length(file) != 1L || !nzchar(file)) stop("`file` must be a single valid file path.", call. = FALSE)
if (is.character(append) && length(append) == 1L && file.exists(append)) {
old <- paste(readLines(append, warn = FALSE, encoding = "UTF-8"), collapse = "\n")
if (grepl("</body>", old, fixed = TRUE)) {
document <- sub("</body>", paste0("<div class=\"table-separator\"></div>", table_block, "</body>"), old, fixed = TRUE)
file <- append
}
}
writeLines(enc2utf8(document), file, useBytes = TRUE)
# Build a flat data frame for Word/Excel directly from the already
# calculated display cells. This is local by design; no separate tabledata()
# helper or HTML parsing is required.
clean_export_text <- function(x) {
x <- as.character(x)
x <- gsub("<br\\s*/?>", " ", x, ignore.case = TRUE)
x <- gsub("<[^>]+>", "", x)
x <- gsub("<", "<", x, fixed = TRUE)
x <- gsub(">", ">", x, fixed = TRUE)
x <- gsub(""", "\"", x, fixed = TRUE)
x <- gsub("'", "'", x, fixed = TRUE)
x <- gsub(" ", " ", x, fixed = TRUE)
x <- gsub("&", "&", x, fixed = TRUE)
x <- gsub("[\r\n\t]+", " ", x)
x <- gsub("\\s+", " ", x)
trimws(x)
}
export_method_titles <- unname(method_labels[methods])
export_method_titles[is.na(export_method_titles) | !nzchar(export_method_titles)] <-
methods[is.na(export_method_titles) | !nzchar(export_method_titles)]
export_column_names <- "Characteristic"
for (i in seq_along(methods)) {
export_column_names <- c(export_column_names, paste0(export_method_titles[i], " ", effect_type, " (95% CI)"))
if (include_p_column) export_column_names <- c(export_column_names, paste0(export_method_titles[i], " p-value"))
}
export_rows <- lapply(rows, function(current) {
characteristic <- if (current$type %in% c("header", "continuous")) {
clean_export_text(current$label)
} else {
paste0(" ", clean_export_text(current$item))
}
output_row <- characteristic
for (i in seq_along(methods)) {
method <- methods[i]
output_row <- c(output_row, clean_export_text(cell_effect_html(method, current)))
if (include_p_column) output_row <- c(output_row, clean_export_text(cell_p_html(method, current)))
}
output_row
})
export_diag_labels <- c("N", "Events", "Variables", "Parameters", "AIC", "BIC",
"McFadden R2", "Nagelkerke R2", "Hosmer-Lemeshow p",
"Pearson GOF p", "AUC", "VIF range")
export_diag_keys <- c("N", "Events", "Variables", "Parameters", "AIC", "BIC",
"McFadden", "Nagelkerke", "Hosmer", "Pearson", "AUC", "VIF")
for (j in seq_along(export_diag_labels)) {
output_row <- export_diag_labels[j]
key <- export_diag_keys[j]
for (i in seq_along(methods)) {
d <- diagnostics[[methods[i]]]
value <- switch(
key,
N = formatC(d["N"], format = "f", digits = 0),
Events = formatC(d["Events"], format = "f", digits = 0),
Variables = formatC(d["Variables"], format = "f", digits = 0),
Parameters = formatC(d["Parameters"], format = "f", digits = 0),
AIC = format_number(d["AIC"], 2),
BIC = format_number(d["BIC"], 2),
McFadden = format_number(d["McFadden"], 3),
Nagelkerke = format_number(d["Nagelkerke"], 3),
Hosmer = format_p_plain(d["Hosmer"]),
Pearson = format_p_plain(d["Pearson"]),
AUC = format_number(d["AUC"], 3),
VIF = if (all(is.finite(d[c("VIF_min", "VIF_max")]))) {
paste0(format_number(d["VIF_min"], 2), " - ", format_number(d["VIF_max"], 2))
} else ""
)
output_row <- c(output_row, clean_export_text(value))
if (include_p_column) output_row <- c(output_row, "")
}
export_rows[[length(export_rows) + 1L]] <- output_row
}
if (length(export_rows)) {
table_df <- as.data.frame(do.call(rbind, export_rows), stringsAsFactors = FALSE, check.names = FALSE)
} else {
table_df <- as.data.frame(matrix(character(), nrow = 0L, ncol = length(export_column_names)),
stringsAsFactors = FALSE, check.names = FALSE)
}
names(table_df) <- make.unique(export_column_names, sep = "_")
rownames(table_df) <- NULL
output <- list(
rows = rows,
data = table_df,
metadata = vars,
by = by_name,
event = event_level,
effect = effect_type,
methods = methods,
global = global,
rvrow = reverse_rows,
selected = selected_by_method,
models = fitted_by_method,
diagnostics = diagnostics,
selection_details = selection_details,
raw = if (isTRUE(raw)) extracted else NULL,
html = document,
table_html = table_block,
blocks = blocks,
file = normalizePath(file, winslash = "/", mustWork = TRUE),
call = match.call()
)
class(output) <- "r4vn_tabmulti"
if (isTRUE(show)) {
viewer <- getOption("viewer")
if (is.function(viewer)) viewer(output$file) else utils::browseURL(output$file)
}
invisible(output)
}
#' Print or Reopen an R4VN Model-Comparison Table
#'
#' Opens the HTML file stored in an \code{r4vn_tabmulti} object.
#'
#' @param x An object created by \code{tabmulti()}.
#' @param ... Additional arguments currently ignored.
#' @return The input object, invisibly.
#' @keywords internal
#' @method print r4vn_tabmulti
#' @export
print.r4vn_tabmulti <- function(x, ...) {
viewer <- getOption("viewer")
if (is.function(viewer)) viewer(x$file) else utils::browseURL(x$file)
invisible(x)
}
# Version marker used to verify that the current source file has been loaded.
attr(tabmulti, "r4vn_version") <- "tabmulti-1.1.0-2026-08-03"
Try the R4VN package in your browser
Any scripts or data that you put into this service are public.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.