Nothing
R/tabmeta.R
In R4VN: Health Data Analysis and Publication-Ready Reporting
Defines functions
.r4vn_meta_subgroups
.r4vn_meta_metric_ci_text
.r4vn_meta_heterogeneity_ci
.r4vn_meta_heterogeneity
.r4vn_meta_binary_table
.r4vn_meta_study_table
.r4vn_meta_term_labels
.r4vn_meta_group_values
.r4vn_meta_group_spec
.r4vn_meta_variable_label
.r4vn_meta_weights
.r4vn_meta_prediction
.r4vn_meta_model
.r4vn_meta_make_es
.r4vn_meta_ci_text_compact
.r4vn_meta_fmt_compact
.r4vn_meta_ci_text
.r4vn_meta_p
.r4vn_meta_fmt
.r4vn_meta_back_fun
.r4vn_meta_back
.r4vn_meta_transform
.r4vn_meta_config
.r4vn_meta_effect_type
.r4vn_meta_test_method
.r4vn_meta_profile
.r4vn_meta_varnames
.r4vn_meta_reg_names
.r4vn_meta_varname
.r4vn_meta_eval
.r4vn_meta_escape
.r4vn_meta_num1
.r4vn_meta_flag
.r4vn_meta_need
# ============================================================================
# R4VN tabmeta: publication-ready meta-analysis
# Backend: metafor (optional Suggests dependency)
# ============================================================================
.r4vn_meta_need <- function(package = "metafor") {
if (!requireNamespace(package, quietly = TRUE)) {
stop(
"Package `", package, "` is required for this analysis. Install it with ",
"install.packages(\"", package, "\").",
call. = FALSE
)
}
invisible(TRUE)
}
.r4vn_meta_flag <- function(x, name) {
if (!is.logical(x) || length(x) != 1L || is.na(x)) {
stop("`", name, "` must be TRUE or FALSE.", call. = FALSE)
}
x
}
.r4vn_meta_num1 <- function(x, name, lower = -Inf, upper = Inf,
inclusive_lower = TRUE, inclusive_upper = TRUE) {
if (!is.numeric(x) || length(x) != 1L || is.na(x) || !is.finite(x)) {
stop("`", name, "` must be one finite numeric value.", call. = FALSE)
}
ok_lower <- if (inclusive_lower) x >= lower else x > lower
ok_upper <- if (inclusive_upper) x <= upper else x < upper
if (!ok_lower || !ok_upper) {
stop("`", name, "` is outside the allowed range.", call. = FALSE)
}
x
}
.r4vn_meta_escape <- function(x) {
x <- as.character(x)
x <- gsub("&", "&", x, fixed = TRUE)
x <- gsub("<", "<", x, fixed = TRUE)
x <- gsub(">", ">", x, fixed = TRUE)
x <- gsub('"', """, x, fixed = TRUE)
gsub("'", "'", x, fixed = TRUE)
}
.r4vn_meta_eval <- function(expr, data, env, arg, required = FALSE) {
if (identical(expr, quote(NULL)) || is.null(expr)) {
if (required) stop("`", arg, "` is required.", call. = FALSE)
return(NULL)
}
value <- tryCatch(
eval(expr, envir = data, enclos = env),
error = function(e) tryCatch(eval(expr, envir = env),
error = function(e2) NULL)
)
if (is.null(value) && required) {
stop("Could not evaluate `", arg, "`.", call. = FALSE)
}
if (!is.null(value) && length(value) == 1L && nrow(data) > 1L &&
is.character(value) && value %in% names(data)) {
value <- data[[value]]
}
value
}
.r4vn_meta_varname <- function(expr, data, env, arg, allow_null = TRUE) {
if (identical(expr, quote(NULL)) || is.null(expr)) {
if (allow_null) return(NULL)
stop("`", arg, "` is required.", call. = FALSE)
}
if (is.symbol(expr)) {
nm <- as.character(expr)
if (nm %in% names(data)) return(nm)
}
value <- tryCatch(eval(expr, envir = data, enclos = env),
error = function(e) NULL)
if (is.character(value) && length(value) == 1L && value %in% names(data)) {
return(value)
}
stop("`", arg, "` must identify one variable in `data`.", call. = FALSE)
}
.r4vn_meta_reg_names <- function(reg, data) {
if (is.null(reg)) return(character())
if (inherits(reg, "r4vn_vars")) {
out <- as.character(reg$variable)
} else if (is.character(reg)) {
out <- as.character(reg)
} else {
stop("`reg` must be created using `vars()` or be a character vector.",
call. = FALSE)
}
out <- unique(out[nzchar(out)])
absent <- setdiff(out, names(data))
if (length(absent)) {
stop("Variables in `reg` not found in `data`: ",
paste(absent, collapse = ", "), ".", call. = FALSE)
}
out
}
.r4vn_meta_varnames <- function(expr, data, env, arg, allow_null = TRUE) {
if (identical(expr, quote(NULL)) || is.null(expr)) {
if (allow_null) return(character())
stop("`", arg, "` is required.", call. = FALSE)
}
if (is.symbol(expr)) {
nm <- as.character(expr)
if (nm %in% names(data)) return(nm)
}
value <- tryCatch(eval(expr, envir = data, enclos = env),
error = function(e) NULL)
if (inherits(value, "r4vn_vars")) {
out <- as.character(value$variable)
} else if (is.character(value)) {
out <- as.character(value)
} else {
stop("`", arg, "` must identify variables in `data` or use `vars()`.",
call. = FALSE)
}
out <- unique(out[!is.na(out) & nzchar(out)])
absent <- setdiff(out, names(data))
if (length(absent)) {
stop("Variables in `", arg, "` not found in `data`: ",
paste(absent, collapse = ", "), ".", call. = FALSE)
}
out
}
.r4vn_meta_profile <- function(profile, full, prediction, bias, leaveout,
influence, plot) {
profile <- match.arg(profile, c("auto", "brief", "full", "custom"))
if (!is.null(full)) {
.r4vn_meta_flag(full, "full")
profile <- if (isTRUE(full)) "full" else "custom"
}
defaults <- switch(
profile,
auto = list(prediction = TRUE, bias = TRUE, leaveout = TRUE,
influence = TRUE, plot = TRUE),
brief = list(prediction = FALSE, bias = FALSE, leaveout = FALSE,
influence = FALSE, plot = TRUE),
full = list(prediction = TRUE, bias = TRUE, leaveout = TRUE,
influence = TRUE, plot = TRUE),
custom = list(prediction = FALSE, bias = FALSE, leaveout = FALSE,
influence = FALSE, plot = FALSE)
)
supplied <- list(prediction = prediction, bias = bias, leaveout = leaveout,
influence = influence, plot = plot)
for (nm in names(defaults)) {
if (!is.null(supplied[[nm]])) {
.r4vn_meta_flag(supplied[[nm]], nm)
defaults[[nm]] <- supplied[[nm]]
}
}
c(defaults, list(profile = profile, full = identical(profile, "full")))
}
.r4vn_meta_test_method <- function(small, k, random = TRUE) {
small <- match.arg(tolower(as.character(small)[1L]),
c("auto", "adhoc", "knha", "t", "z"))
if (!isTRUE(random)) return("z")
if (identical(small, "auto")) {
if (k <= 10L) "adhoc" else "z"
} else small
}
.r4vn_meta_effect_type <- function(or, rr, rd, hr, irr, md, smd, prop,
rate, cor, effect_present,
inferred = NULL) {
flags <- c(OR = or, RR = rr, RD = rd, HR = hr, IRR = irr, MD = md,
SMD = smd, PROP = prop, RATE = rate, COR = cor)
selected <- names(flags)[flags]
if (length(selected) > 1L) {
stop("Choose only one effect option: `or`, `rr`, `rd`, `hr`, `irr`, ",
"`md`, `smd`, `prop`, `rate`, or `cor`.", call. = FALSE)
}
if (!length(selected)) {
if (effect_present) return("GENERIC")
if (!is.null(inferred) && length(inferred)) return(inferred[1L])
stop("Choose an effect type such as `or=TRUE`, `rr=TRUE`, `hr=TRUE`, ",
"`md=TRUE`, `smd=TRUE`, `prop=TRUE`, `rate=TRUE`, or `cor=TRUE`.",
call. = FALSE)
}
selected
}
.r4vn_meta_config <- function(type, transform = NULL) {
if (type %in% c("OR", "RR", "HR", "IRR")) {
return(list(measure = type, ref = 1, label = type, scale = "log"))
}
if (type == "RD") {
return(list(measure = "RD", ref = 0, label = "Risk difference",
scale = "identity"))
}
if (type == "MD") {
return(list(measure = "MD", ref = 0, label = "Mean difference",
scale = "identity"))
}
if (type == "SMD") {
return(list(measure = "SMD", ref = 0,
label = "Standardized mean difference", scale = "identity"))
}
if (type == "COR") {
return(list(measure = "ZCOR", ref = 0, label = "Correlation",
scale = "zcor"))
}
if (type == "RATE") {
return(list(measure = "IRLN", ref = 0, label = "Incidence rate",
scale = "log"))
}
if (type == "PROP") {
if (is.null(transform)) transform <- "logit"
transform <- tolower(as.character(transform)[1L])
transform <- switch(
transform,
logit = "logit",
arcsine = "arcsine",
asin = "arcsine",
ft = "ft",
`freeman-tukey` = "ft",
none = "none",
stop("`transform` for proportions must be `logit`, `arcsine`, `ft`, or `none`.",
call. = FALSE)
)
measure <- switch(transform, logit = "PLO", arcsine = "PAS",
ft = "PFT", none = "PR")
return(list(measure = measure, ref = 0, label = "Proportion",
scale = transform))
}
list(measure = "GEN", ref = 0, label = "Effect", scale = "identity")
}
.r4vn_meta_transform <- function(x, config) {
if (config$scale == "log") {
if (any(x <= 0, na.rm = TRUE)) {
stop("Ratio/rate estimates must be greater than 0.", call. = FALSE)
}
return(log(x))
}
if (config$scale == "zcor") {
if (any(abs(x) >= 1, na.rm = TRUE)) {
stop("Correlation estimates must be strictly between -1 and 1.",
call. = FALSE)
}
return(atanh(x))
}
x
}
.r4vn_meta_back <- function(x, config, ni = NULL) {
if (config$scale == "log") return(exp(x))
if (config$scale == "zcor") return(tanh(x))
if (config$scale == "logit") return(stats::plogis(x))
if (config$scale == "arcsine") return(sin(x)^2)
if (config$scale == "ft") {
.r4vn_meta_need()
if (is.null(ni)) return(rep(NA_real_, length(x)))
if (length(x) == length(ni) && length(x) > 1L) {
return(mapply(
function(z, nn) metafor::transf.ipft(z, targs = list(ni = nn)),
x, ni, USE.NAMES = FALSE
))
}
return(metafor::transf.ipft.hm(x, targs = list(ni = ni)))
}
x
}
.r4vn_meta_back_fun <- function(x) {
config <- x$config
ni <- x$analysis_data$.n_for_back
if (config$scale == "log") return(exp)
if (config$scale == "zcor") return(tanh)
if (config$scale == "logit") return(stats::plogis)
if (config$scale == "arcsine") return(function(z) sin(z)^2)
if (config$scale == "ft") {
return(function(z) metafor::transf.ipft.hm(z, targs = list(ni = ni)))
}
identity
}
.r4vn_meta_fmt <- function(x, digits = 2) {
if (!length(x)) return(character(0))
out <- rep("", length(x))
ok <- !is.na(x) & is.finite(x)
if (any(ok)) {
out[ok] <- formatC(
x[ok], format = "f", digits = digits, big.mark = ","
)
}
out
}
.r4vn_meta_p <- function(x, digits = 3) {
if (!length(x) || is.na(x) || !is.finite(x)) return("")
lim <- 10^(-digits)
if (x < lim) paste0("<", formatC(lim, format = "f", digits = digits)) else
formatC(x, format = "f", digits = digits)
}
.r4vn_meta_ci_text <- function(est, lo, hi, digits = 2) {
if (any(!is.finite(c(est, lo, hi)))) return("")
paste0(.r4vn_meta_fmt(est, digits), " (",
.r4vn_meta_fmt(lo, digits), "-", .r4vn_meta_fmt(hi, digits), ")")
}
.r4vn_meta_fmt_compact <- function(x, digits = 2) {
if (!length(x)) return(character())
out <- .r4vn_meta_fmt(x, digits)
finite <- is.finite(x)
scientific <- finite & x != 0 &
(abs(x) >= 1e6 | abs(x) < 10^(-digits))
if (any(scientific)) {
out[scientific] <- formatC(
x[scientific], format = "e", digits = digits
)
}
out
}
.r4vn_meta_ci_text_compact <- function(est, lo, hi, digits = 2) {
if (any(!is.finite(c(est, lo, hi)))) return("")
paste0(
.r4vn_meta_fmt_compact(est, digits), " (",
.r4vn_meta_fmt_compact(lo, digits), "-",
.r4vn_meta_fmt_compact(hi, digits), ")"
)
}
.r4vn_meta_make_es <- function(type, config, values, cc, zero) {
.r4vn_meta_need()
drop00 <- identical(zero, "exclude")
if (type %in% c("OR", "RR", "RD")) {
ai <- values$event1
ci0 <- values$event0
n1 <- values$n1
n0 <- values$n0
if (any(ai < 0 | ci0 < 0 | n1 <= 0 | n0 <= 0, na.rm = TRUE) ||
any(ai > n1 | ci0 > n0, na.rm = TRUE)) {
stop("Binary counts must satisfy 0 <= event <= n and n > 0.",
call. = FALSE)
}
return(metafor::escalc(
measure = config$measure,
ai = ai, bi = n1 - ai,
ci = ci0, di = n0 - ci0,
add = cc, to = if (cc == 0) "none" else "only0",
drop00 = drop00
))
}
if (type %in% c("MD", "SMD")) {
if (any(values$n1 <= 1 | values$n0 <= 1, na.rm = TRUE) ||
any(values$sd1 < 0 | values$sd0 < 0, na.rm = TRUE)) {
stop("Continuous input requires n1,n0 > 1 and non-negative SDs.",
call. = FALSE)
}
return(metafor::escalc(
measure = config$measure,
m1i = values$mean1, sd1i = values$sd1, n1i = values$n1,
m2i = values$mean0, sd2i = values$sd0, n2i = values$n0
))
}
if (type == "PROP") {
if (any(values$event < 0 | values$n <= 0 | values$event > values$n,
na.rm = TRUE)) {
stop("Proportion input requires 0 <= event <= n and n > 0.",
call. = FALSE)
}
return(metafor::escalc(
measure = config$measure, xi = values$event, ni = values$n,
add = cc, to = if (cc == 0) "none" else "only0"
))
}
if (type == "RATE") {
if (any(values$event < 0 | values$time <= 0, na.rm = TRUE)) {
stop("Rate input requires event >= 0 and person-time > 0.",
call. = FALSE)
}
return(metafor::escalc(
measure = "IRLN", xi = values$event, ti = values$time,
add = cc, to = if (cc == 0) "none" else "only0"
))
}
if (type == "IRR") {
if (any(values$event1 < 0 | values$event0 < 0 |
values$time1 <= 0 | values$time0 <= 0, na.rm = TRUE)) {
stop("IRR input requires non-negative events and positive person-time.",
call. = FALSE)
}
return(metafor::escalc(
measure = "IRR",
x1i = values$event1, x2i = values$event0,
t1i = values$time1, t2i = values$time0,
add = cc, to = if (cc == 0) "none" else "only0"
))
}
if (type == "COR") {
if (any(abs(values$effect) >= 1, na.rm = TRUE) ||
any(values$n <= 3, na.rm = TRUE)) {
stop("Correlation input requires -1 < effect < 1 and n > 3.",
call. = FALSE)
}
return(metafor::escalc(
measure = "ZCOR", ri = values$effect, ni = values$n
))
}
NULL
}
.r4vn_meta_model <- function(yi, vi, method, hk = NULL, ci, fixed = FALSE,
mods = NULL, test = NULL) {
.r4vn_meta_need()
if (is.null(test)) {
test <- if (!fixed && isTRUE(hk)) "knha" else "z"
}
test <- .r4vn_meta_test_method(test, length(yi), random = !fixed)
args <- list(
yi = yi,
vi = vi,
method = if (fixed) "FE" else method,
test = test,
level = ci * 100
)
if (!is.null(mods)) args$mods <- mods
do.call(metafor::rma.uni, args)
}
.r4vn_meta_prediction <- function(model, config, ci, ni = NULL) {
pred <- stats::predict(model, level = ci * 100)
list(
estimate = .r4vn_meta_back(as.numeric(pred$pred), config, ni),
lower = .r4vn_meta_back(as.numeric(pred$ci.lb), config, ni),
upper = .r4vn_meta_back(as.numeric(pred$ci.ub), config, ni),
pi_lower = if (!is.null(pred$pi.lb))
.r4vn_meta_back(as.numeric(pred$pi.lb), config, ni) else NA_real_,
pi_upper = if (!is.null(pred$pi.ub))
.r4vn_meta_back(as.numeric(pred$pi.ub), config, ni) else NA_real_
)
}
.r4vn_meta_weights <- function(model) {
w <- tryCatch(stats::weights(model), error = function(e) NULL)
if (is.null(w)) return(rep(NA_real_, model$k))
as.numeric(w)
}
.r4vn_meta_variable_label <- function(x, fallback) {
label <- attr(x, "label", exact = TRUE)
if (is.null(label) || !length(label) || is.na(label[1L]) ||
!nzchar(trimws(as.character(label)[1L]))) {
return(as.character(fallback)[1L])
}
as.character(label)[1L]
}
.r4vn_meta_group_spec <- function(x) {
observed <- x[!is.na(x)]
raw <- if (is.factor(x)) {
levels(droplevels(x))
} else {
unique(as.character(observed))
}
raw <- raw[raw %in% as.character(observed)]
display <- raw
value_labels <- attr(x, "labels", exact = TRUE)
if (is.null(value_labels)) {
value_labels <- attr(x, "value.labels", exact = TRUE)
}
if (!is.null(value_labels) && length(value_labels) &&
!is.null(names(value_labels))) {
labelled_values <- as.character(unname(value_labels))
matched <- match(raw, labelled_values)
use <- !is.na(matched) & nzchar(names(value_labels)[matched])
display[use] <- names(value_labels)[matched[use]]
}
duplicate <- duplicated(display) | duplicated(display, fromLast = TRUE)
display[duplicate] <- paste0(display[duplicate], " (", raw[duplicate], ")")
list(raw = raw, display = display)
}
.r4vn_meta_group_values <- function(x) {
spec <- .r4vn_meta_group_spec(x)
out <- as.character(x)
matched <- match(out, spec$raw)
use <- !is.na(matched)
out[use] <- spec$display[matched[use]]
out[is.na(x)] <- NA_character_
out
}
.r4vn_meta_term_labels <- function(terms, data, variables) {
out <- as.character(terms)
normalized <- tolower(gsub("[()[:space:]]", "", out))
out[normalized %in% c("intrcpt", "intercept")] <- "Intercept"
variables <- unique(as.character(variables))
variables <- variables[variables %in% names(data)]
if (!length(variables)) return(out)
variables <- variables[order(nchar(variables), decreasing = TRUE)]
for (nm in variables) {
variable_label <- .r4vn_meta_variable_label(data[[nm]], nm)
exact <- terms == nm
out[exact] <- variable_label
prefixed <- !exact & startsWith(as.character(terms), nm)
if (any(prefixed)) {
suffix <- substring(as.character(terms)[prefixed], nchar(nm) + 1L)
suffix <- sub("^[.:_]+", "", suffix)
out[prefixed] <- ifelse(
nzchar(suffix), paste0(variable_label, ": ", suffix), variable_label
)
}
}
out
}
.r4vn_meta_study_table <- function(d, model, config, digits,
subgroup = NULL, subgroup_title = "Subgroup",
ci = 0.95) {
back <- .r4vn_meta_back_fun(list(config = config, analysis_data = d))
zcrit <- stats::qnorm(1 - (1 - ci) / 2)
est <- back(d$yi)
lo <- back(d$yi - zcrit * sqrt(d$vi))
hi <- back(d$yi + zcrit * sqrt(d$vi))
wt <- .r4vn_meta_weights(model)
if (length(wt) != nrow(d)) wt <- rep(NA_real_, nrow(d))
out <- data.frame(
Study = d$.study,
stringsAsFactors = FALSE, check.names = FALSE
)
if (all(c(".event1", ".n1", ".event0", ".n0") %in% names(d))) {
out[["Group 1"]] <- paste0(
.r4vn_meta_fmt(d$.event1, 0), "/", .r4vn_meta_fmt(d$.n1, 0)
)
out[["Group 0"]] <- paste0(
.r4vn_meta_fmt(d$.event0, 0), "/", .r4vn_meta_fmt(d$.n0, 0)
)
} else if (all(c(".event1", ".time1", ".event0", ".time0") %in% names(d))) {
out[["Group 1"]] <- paste0(
.r4vn_meta_fmt(d$.event1, 0), "/",
.r4vn_meta_fmt(d$.time1, digits)
)
out[["Group 0"]] <- paste0(
.r4vn_meta_fmt(d$.event0, 0), "/",
.r4vn_meta_fmt(d$.time0, digits)
)
} else if (all(c(".mean1", ".sd1", ".n1", ".mean0", ".sd0", ".n0") %in%
names(d))) {
out[["Group 1"]] <- paste0(
.r4vn_meta_fmt(d$.mean1, digits), " (",
.r4vn_meta_fmt(d$.sd1, digits), "); n=",
.r4vn_meta_fmt(d$.n1, 0)
)
out[["Group 0"]] <- paste0(
.r4vn_meta_fmt(d$.mean0, digits), " (",
.r4vn_meta_fmt(d$.sd0, digits), "); n=",
.r4vn_meta_fmt(d$.n0, 0)
)
} else if (all(c(".event", ".n") %in% names(d))) {
out[["Events/Total"]] <- paste0(
.r4vn_meta_fmt(d$.event, 0), "/", .r4vn_meta_fmt(d$.n, 0)
)
} else if (all(c(".event", ".time") %in% names(d))) {
out[["Events/Person-time"]] <- paste0(
.r4vn_meta_fmt(d$.event, 0), "/",
.r4vn_meta_fmt(d$.time, digits)
)
}
out[["Effect (95% CI)"]] <- mapply(
.r4vn_meta_ci_text, est, lo, hi,
MoreArgs = list(digits = digits), USE.NAMES = FALSE
)
out[["Weight"]] <- ifelse(
is.finite(wt), paste0(.r4vn_meta_fmt(wt, 1), "%"), ""
)
if (!is.null(subgroup)) {
subgroup_data <- data.frame(subgroup, stringsAsFactors = FALSE,
check.names = FALSE)
names(subgroup_data) <- as.character(subgroup_title)[1L]
out <- cbind(subgroup_data, out)
}
out
}
.r4vn_meta_binary_table <- function(d) {
required <- c(".event1", ".n1", ".event0", ".n0")
if (!all(required %in% names(d))) return(NULL)
a <- d$.event1
b <- d$.n1 - d$.event1
c0 <- d$.event0
d0 <- d$.n0 - d$.event0
data.frame(
Study = d$.study,
`a: event, group 1` = a,
`b: non-event, group 1` = b,
`c: event, group 0` = c0,
`d: non-event, group 0` = d0,
`Group 1 total` = d$.n1,
`Group 0 total` = d$.n0,
stringsAsFactors = FALSE, check.names = FALSE
)
}
.r4vn_meta_heterogeneity <- function(model, p_digits = 3) {
data.frame(
Statistic = c("Studies", "Q", "df", "p heterogeneity", "I-squared",
"H-squared", "tau-squared", "tau"),
Value = c(
as.character(model$k),
.r4vn_meta_fmt(model$QE, 2),
as.character(model$k - model$p),
.r4vn_meta_p(model$QEp, p_digits),
paste0(.r4vn_meta_fmt(model$I2, 1), "%"),
.r4vn_meta_fmt(model$H2, 2),
.r4vn_meta_fmt(model$tau2, 4),
.r4vn_meta_fmt(sqrt(model$tau2), 4)
),
stringsAsFactors = FALSE, check.names = FALSE
)
}
.r4vn_meta_heterogeneity_ci <- function(model, ci = 0.95) {
empty <- list(
i2 = c(estimate = NA_real_, lower = NA_real_, upper = NA_real_),
tau2 = c(estimate = NA_real_, lower = NA_real_, upper = NA_real_)
)
result <- tryCatch(
suppressWarnings(stats::confint(model, level = ci * 100)),
error = function(e) NULL
)
if (is.null(result)) return(empty)
random <- if (is.list(result) && !is.null(result$random)) {
result$random
} else if (is.matrix(result) || is.data.frame(result)) {
result
} else NULL
if (is.null(random)) return(empty)
random <- as.data.frame(random, check.names = FALSE)
if (!nrow(random) || ncol(random) < 3L) return(empty)
row_key <- tolower(gsub("[^a-zA-Z0-9]", "", rownames(random)))
column_key <- tolower(gsub("[^a-zA-Z0-9]", "", names(random)))
estimate_column <- match("estimate", column_key)
lower_column <- match("cilb", column_key)
upper_column <- match("ciub", column_key)
if (anyNA(c(estimate_column, lower_column, upper_column))) {
estimate_column <- 1L
lower_column <- 2L
upper_column <- 3L
}
extract <- function(key) {
row <- match(key, row_key)
if (is.na(row)) return(c(estimate = NA_real_, lower = NA_real_,
upper = NA_real_))
values <- suppressWarnings(as.numeric(random[
row, c(estimate_column, lower_column, upper_column), drop = TRUE
]))
stats::setNames(values, c("estimate", "lower", "upper"))
}
list(i2 = extract("i2"), tau2 = extract("tau2"))
}
.r4vn_meta_metric_ci_text <- function(estimate, interval, digits = 2,
suffix = "") {
if (!is.finite(estimate)) return("")
if (length(interval) >= 3L &&
all(is.finite(interval[c("lower", "upper")]))) {
return(paste0(
.r4vn_meta_fmt(estimate, digits), suffix, " (",
.r4vn_meta_fmt(interval[["lower"]], digits), suffix, "-",
.r4vn_meta_fmt(interval[["upper"]], digits), suffix, ")"
))
}
paste0(.r4vn_meta_fmt(estimate, digits), suffix)
}
.r4vn_meta_subgroups <- function(d, by_name, method, test, ci, config,
digits, p_digits, fixed = FALSE,
include_interpretation = FALSE) {
if (is.null(by_name)) return(NULL)
g <- d[[by_name]]
group_spec <- .r4vn_meta_group_spec(g)
lev <- group_spec$raw
display_lev <- group_spec$display
variable_label <- .r4vn_meta_variable_label(g, by_name)
if (length(lev) < 2L) return(NULL)
fits <- list()
indices <- list()
rows <- list()
retained_raw <- retained_display <- character()
for (i in seq_along(lev)) {
z <- lev[i]
display_z <- display_lev[i]
idx <- !is.na(g) & as.character(g) == z
if (sum(idx) < 2L) next
fit <- .r4vn_meta_model(d$yi[idx], d$vi[idx], method, ci = ci,
test = test, fixed = fixed)
fits[[display_z]] <- fit
indices[[display_z]] <- idx
retained_raw <- c(retained_raw, z)
retained_display <- c(retained_display, display_z)
pr <- .r4vn_meta_prediction(fit, config, ci, d$.n_for_back[idx])
heterogeneity_ci <- .r4vn_meta_heterogeneity_ci(fit, ci)
fit_test <- tolower(as.character(fit$test %||% test)[1L])
test_symbol <- if (fit_test %in% c("knha", "adhoc", "t")) "t" else "z"
effect_statistic <- as.numeric(fit$zval)[1L]
effect_df <- if (test_symbol == "t") {
as.numeric((fit$ddf %||% (fit$k - fit$p))[1L])
} else NA_real_
effect_p <- as.numeric(fit$pval)[1L]
prediction_interval <- if (!fixed && is.finite(pr$pi_lower) &&
is.finite(pr$pi_upper)) {
paste0(
.r4vn_meta_fmt(pr$pi_lower, digits), "-",
.r4vn_meta_fmt(pr$pi_upper, digits)
)
} else ""
subgroup_row <- data.frame(
Subgroup = display_z,
Studies = fit$k,
`Effect (95% CI)` = .r4vn_meta_ci_text(
pr$estimate, pr$lower, pr$upper, digits
),
`Test statistic` = if (is.finite(effect_statistic)) {
paste0(test_symbol, " = ", .r4vn_meta_fmt(effect_statistic, 2))
} else "",
`df effect` = if (is.finite(effect_df)) {
.r4vn_meta_fmt(effect_df, 0)
} else "",
`p-value effect` = .r4vn_meta_p(effect_p, p_digits),
`Prediction interval` = prediction_interval,
`Q heterogeneity` = .r4vn_meta_fmt(fit$QE, 2),
`df heterogeneity` = as.character(max(0, fit$k - fit$p)),
`p-value heterogeneity` = .r4vn_meta_p(fit$QEp, p_digits),
`I-squared (95% CI)` = .r4vn_meta_metric_ci_text(
fit$I2, heterogeneity_ci$i2, digits = 1, suffix = "%"
),
`Tau-squared (95% CI)` = .r4vn_meta_metric_ci_text(
fit$tau2, heterogeneity_ci$tau2, digits = 4
),
stringsAsFactors = FALSE, check.names = FALSE
)
confidence_percent <- round(ci * 100)
names(subgroup_row)[names(subgroup_row) == "Effect (95% CI)"] <-
paste0("Effect (", confidence_percent, "% CI)")
names(subgroup_row)[names(subgroup_row) == "Prediction interval"] <-
paste0(confidence_percent, "% prediction interval")
names(subgroup_row)[names(subgroup_row) == "I-squared (95% CI)"] <-
paste0("I-squared (", confidence_percent, "% CI)")
names(subgroup_row)[names(subgroup_row) == "Tau-squared (95% CI)"] <-
paste0("Tau-squared (", confidence_percent, "% CI)")
if (isTRUE(include_interpretation)) {
effect_significant <- is.finite(effect_p) && effect_p < (1 - ci)
heterogeneity_significant <- is.finite(fit$QEp) &&
fit$QEp < (1 - ci)
effect_conclusion <- if (!is.finite(effect_p)) {
"the pooled-effect significance test was unavailable"
} else if (effect_significant) {
"the pooled effect was statistically significant"
} else {
"the pooled effect was not statistically significant"
}
heterogeneity_conclusion <- if (!is.finite(fit$QEp)) {
"the heterogeneity significance test was unavailable."
} else if (heterogeneity_significant) {
"heterogeneity was statistically significant."
} else {
"heterogeneity was not statistically significant."
}
subgroup_row[["Statistically significant"]] <- if (
is.finite(effect_p)
) if (effect_significant) "Yes" else "No" else ""
subgroup_row[["Conclusion"]] <- paste0(
toupper(substr(effect_conclusion, 1L, 1L)),
substring(effect_conclusion, 2L), "; ",
heterogeneity_conclusion
)
}
rows[[display_z]] <- subgroup_row
}
table <- if (length(rows)) do.call(rbind, rows) else NULL
if (!is.null(table)) rownames(table) <- NULL
complete_group <- !is.na(g)
moddata <- data.frame(
.g = factor(as.character(g[complete_group]))
)
mm <- stats::model.matrix(~ .g, data = moddata)
testfit <- if (ncol(mm) > 1L) tryCatch(
.r4vn_meta_model(
d$yi[complete_group], d$vi[complete_group], method, ci = ci,
mods = mm[, -1L, drop = FALSE], test = test, fixed = fixed
),
error = function(e) NULL
) else NULL
test <- if (is.null(testfit)) NULL else data.frame(
Statistic = c("Q between subgroups", "df", "p"),
Value = c(
.r4vn_meta_fmt(testfit$QM, 2),
as.character(testfit$m),
.r4vn_meta_p(testfit$QMp, p_digits)
),
stringsAsFactors = FALSE, check.names = FALSE
)
list(
fits = fits, indices = indices, table = table, test = test,
test_model = testfit, levels = retained_display,
raw_levels = retained_raw, variable = by_name, label = variable_label
)
}
.r4vn_meta_subgroup_transpose <- function(table) {
if (is.null(table) || !is.data.frame(table) || !nrow(table) ||
!"Subgroup" %in% names(table)) return(NULL)
metrics <- setdiff(names(table), "Subgroup")
out <- data.frame(
Statistic = metrics, stringsAsFactors = FALSE, check.names = FALSE
)
for (i in seq_len(nrow(table))) {
group_name <- as.character(table$Subgroup[i])
values <- vapply(
table[i, metrics, drop = FALSE],
function(value) {
if (!length(value) || is.na(value[1L])) "" else as.character(value[1L])
},
character(1)
)
out[[group_name]] <- unname(values)
}
out
}
.r4vn_meta_significance_tests <- function(model, overall, test_method,
ci = 0.95, p_digits = 3,
include_interpretation = FALSE) {
alpha <- 1 - ci
pooled_stat <- if (!is.null(model$zval) && length(model$zval)) {
as.numeric(model$zval[1L])
} else {
NA_real_
}
pooled_df <- if (test_method %in% c("t", "knha", "adhoc")) {
candidate <- model$ddf %||% (model$k - model$p)
as.numeric(candidate[1L])
} else {
NA_real_
}
pooled_p <- as.numeric(overall$p)[1L]
heterogeneity_p <- as.numeric(model$QEp)[1L]
pooled_sig <- is.finite(pooled_p) && pooled_p < alpha
heterogeneity_sig <- is.finite(heterogeneity_p) &&
heterogeneity_p < alpha
pooled_symbol <- if (test_method %in% c("t", "knha", "adhoc")) {
"t"
} else {
"z"
}
out <- data.frame(
Test = c("Pooled effect", "Heterogeneity (Cochran's Q)"),
`Test statistic` = c(
if (is.finite(pooled_stat)) {
paste0(pooled_symbol, " = ", .r4vn_meta_fmt(pooled_stat, 2))
} else "",
if (is.finite(model$QE)) {
paste0("Q = ", .r4vn_meta_fmt(model$QE, 2))
} else ""
),
df = c(
if (is.finite(pooled_df)) .r4vn_meta_fmt(pooled_df, 0) else "",
.r4vn_meta_fmt(max(0, model$k - model$p), 0)
),
`p-value` = c(
.r4vn_meta_p(pooled_p, p_digits),
.r4vn_meta_p(heterogeneity_p, p_digits)
),
stringsAsFactors = FALSE, check.names = FALSE
)
if (!isTRUE(include_interpretation)) return(out)
out[["Statistically significant"]] <- c(
if (is.finite(pooled_p)) if (pooled_sig) "Yes" else "No" else "",
if (is.finite(heterogeneity_p)) {
if (heterogeneity_sig) "Yes" else "No"
} else ""
)
out[["Conclusion"]] <- c(
if (!is.finite(pooled_p)) {
"Statistical significance of the pooled estimate was not available."
} else if (pooled_sig) {
"The pooled estimate is statistically significant."
} else {
"The pooled estimate is not statistically significant."
},
if (!is.finite(heterogeneity_p)) {
"Statistical significance of heterogeneity was not available."
} else if (heterogeneity_sig) {
"There is statistically significant between-study heterogeneity."
} else {
"There is no statistically significant evidence of between-study heterogeneity."
}
)
out
}
.r4vn_meta_subgroups_all <- function(d, by_names, method, test, ci, config,
digits, p_digits, fixed = FALSE,
include_interpretation = FALSE) {
if (!length(by_names)) return(list(primary = NULL, analyses = list(),
wide = list(), table = NULL,
test = NULL))
analyses <- lapply(by_names, function(nm) {
.r4vn_meta_subgroups(
d, nm, method, test, ci, config, digits, p_digits, fixed = fixed,
include_interpretation = include_interpretation
)
})
names(analyses) <- by_names
analyses <- Filter(Negate(is.null), analyses)
wide <- lapply(analyses, function(result) {
.r4vn_meta_subgroup_transpose(result$table)
})
wide <- Filter(Negate(is.null), wide)
bind_component <- function(component) {
rows <- lapply(names(analyses), function(nm) {
z <- analyses[[nm]][[component]]
if (is.null(z) || !is.data.frame(z) || !nrow(z)) return(NULL)
data.frame(Moderator = analyses[[nm]]$label %||% nm, z,
stringsAsFactors = FALSE,
check.names = FALSE)
})
rows <- Filter(Negate(is.null), rows)
if (!length(rows)) NULL else {
out <- do.call(rbind, rows)
rownames(out) <- NULL
out
}
}
list(
primary = if (length(analyses)) analyses[[1L]] else NULL,
analyses = analyses,
wide = wide,
table = bind_component("table"),
test = bind_component("test")
)
}
.r4vn_meta_regression <- function(d, reg_names, method, test, ci, digits,
p_digits, config = NULL, fixed = FALSE) {
if (!length(reg_names)) return(NULL)
md <- d[, reg_names, drop = FALSE]
keep <- stats::complete.cases(md)
if (sum(keep) < max(3L, length(reg_names) + 2L)) {
return(list(status = "Not enough complete studies for meta-regression."))
}
mm <- stats::model.matrix(~ ., data = md[keep, , drop = FALSE])
if (ncol(mm) < 2L) {
return(list(status = "No estimable moderator terms."))
}
fit <- .r4vn_meta_model(
d$yi[keep], d$vi[keep], method, ci = ci,
mods = mm[, -1L, drop = FALSE], test = test, fixed = fixed
)
sm <- summary(fit)
beta <- as.numeric(sm$beta)
se <- as.numeric(sm$se)
p <- as.numeric(sm$pval)
rn <- rownames(sm$beta)
if (is.null(rn)) rn <- c("Intercept", colnames(mm)[-1L])
display_rn <- .r4vn_meta_term_labels(rn, d, reg_names)
lower_beta <- if (!is.null(sm$ci.lb)) as.numeric(sm$ci.lb) else {
beta - stats::qnorm(1 - (1 - ci) / 2) * se
}
upper_beta <- if (!is.null(sm$ci.ub)) as.numeric(sm$ci.ub) else {
beta + stats::qnorm(1 - (1 - ci) / 2) * se
}
tab <- data.frame(
Moderator = display_rn,
`Beta (95% CI)` = mapply(
.r4vn_meta_ci_text_compact,
beta, lower_beta, upper_beta,
MoreArgs = list(digits = digits), USE.NAMES = FALSE
),
`p-value` = vapply(p, .r4vn_meta_p, character(1), digits = p_digits),
stringsAsFactors = FALSE, check.names = FALSE
)
confidence_percent <- round(ci * 100)
names(tab)[names(tab) == "Beta (95% CI)"] <- paste0(
"Beta (", confidence_percent, "% CI)"
)
if (!is.null(config) && identical(config$scale, "log")) {
tab[["Ratio of effects (95% CI)"]] <- mapply(
.r4vn_meta_ci_text_compact,
exp(beta), exp(lower_beta), exp(upper_beta),
MoreArgs = list(digits = digits), USE.NAMES = FALSE
)
names(tab)[names(tab) == "Ratio of effects (95% CI)"] <- paste0(
"Ratio of effects (", confidence_percent, "% CI)"
)
}
stats_tab <- data.frame(
Statistic = c("QM", "df", "p moderators", "Residual I-squared",
"Residual tau-squared"),
Value = c(
.r4vn_meta_fmt(fit$QM, 2),
as.character(fit$m),
.r4vn_meta_p(fit$QMp, p_digits),
paste0(.r4vn_meta_fmt(fit$I2, 1), "%"),
.r4vn_meta_fmt(fit$tau2, 4)
),
stringsAsFactors = FALSE, check.names = FALSE
)
univariable <- lapply(reg_names, function(nm) {
mdi <- d[, nm, drop = FALSE]
keep_i <- stats::complete.cases(mdi)
if (sum(keep_i) < 3L) return(NULL)
mm_i <- stats::model.matrix(~ ., data = mdi[keep_i, , drop = FALSE])
if (ncol(mm_i) < 2L) return(NULL)
fit_i <- tryCatch(
.r4vn_meta_model(
d$yi[keep_i], d$vi[keep_i], method, ci = ci,
mods = mm_i[, -1L, drop = FALSE], test = test, fixed = fixed
),
error = function(e) NULL
)
if (is.null(fit_i)) return(NULL)
data.frame(
Moderator = .r4vn_meta_variable_label(d[[nm]], nm),
Studies = sum(keep_i),
QM = .r4vn_meta_fmt(fit_i$QM, 2),
df = fit_i$m,
`p-value` = .r4vn_meta_p(fit_i$QMp, p_digits),
`Residual I-squared` = paste0(.r4vn_meta_fmt(fit_i$I2, 1), "%"),
stringsAsFactors = FALSE, check.names = FALSE
)
})
univariable <- Filter(Negate(is.null), univariable)
list(model = fit, table = tab, statistics = stats_tab,
univariable = if (length(univariable)) do.call(rbind, univariable) else NULL,
keep = keep, model_matrix = mm)
}
.r4vn_meta_bias_methods <- function(methods, profile = "auto") {
if (is.null(methods) || !length(methods) || identical(methods, "auto")) {
methods <- if (identical(profile, "full")) {
c("egger", "begg", "trimfill", "failsafe", "selection")
} else c("egger", "begg", "trimfill")
}
methods <- unique(tolower(as.character(methods)))
aliases <- c(`rank` = "begg", `trim-and-fill` = "trimfill",
`fail-safe` = "failsafe", `selmodel` = "selection")
hit <- methods %in% names(aliases)
methods[hit] <- unname(aliases[methods[hit]])
allowed <- c("egger", "begg", "trimfill", "failsafe", "selection")
bad <- setdiff(methods, allowed)
if (length(bad)) {
stop("Unsupported publication-bias method: ", paste(bad, collapse = ", "),
".", call. = FALSE)
}
methods
}
.r4vn_meta_bias <- function(model, k, p_digits, methods, config, ci,
ni = NULL, strict = FALSE) {
rows <- list()
adjusted <- list()
notes <- character()
add_test <- function(test, statistic = "", p = "", result = "Available") {
rows[[length(rows) + 1L]] <<- data.frame(
Test = test, Statistic = statistic, `p-value` = p, Result = result,
stringsAsFactors = FALSE, check.names = FALSE
)
}
safe <- function(expr, label) {
warning_messages <- character()
value <- tryCatch(
withCallingHandlers(
expr,
warning = function(w) {
warning_messages <<- c(warning_messages, conditionMessage(w))
invokeRestart("muffleWarning")
}
),
error = function(e) {
if (isTRUE(strict)) stop(e)
notes <<- c(
notes,
paste0(label, " was unavailable: ", conditionMessage(e))
)
NULL
}
)
if (length(warning_messages)) {
notes <<- c(
notes,
paste0(
label, " warning: ",
paste(unique(warning_messages), collapse = " ")
)
)
}
value
}
egger <- begg <- trimfill <- failsafe <- selection <- NULL
if ("egger" %in% methods) {
if (k < 10L) {
add_test("Egger regression test", result = "Not performed (<10 studies)")
} else {
egger <- safe(metafor::regtest(model, model = "lm", predictor = "sei"),
"Egger test")
if (!is.null(egger)) add_test(
"Egger regression test", .r4vn_meta_fmt(egger$zval, 3),
.r4vn_meta_p(egger$pval, p_digits)
)
}
}
if ("begg" %in% methods) {
if (k < 10L) {
add_test("Begg-Mazumdar rank correlation",
result = "Not performed (<10 studies)")
} else {
begg <- safe(metafor::ranktest(model), "Rank-correlation test")
if (!is.null(begg)) {
stat <- if (!is.null(begg$tau)) begg$tau else begg$zval
add_test("Begg-Mazumdar rank correlation",
.r4vn_meta_fmt(stat, 3),
.r4vn_meta_p(begg$pval, p_digits))
}
}
}
if ("trimfill" %in% methods && k >= 3L) {
trimfill <- safe(metafor::trimfill(model), "Trim-and-fill analysis")
if (!is.null(trimfill)) {
pr <- .r4vn_meta_prediction(trimfill, config, ci, ni)
adjusted[[length(adjusted) + 1L]] <- data.frame(
Method = "Trim-and-fill",
`Adjusted effect (95% CI)` = .r4vn_meta_ci_text(
pr$estimate, pr$lower, pr$upper, 3
),
`Imputed studies` = as.integer(trimfill$k0 %||% 0L),
stringsAsFactors = FALSE, check.names = FALSE
)
}
}
if ("failsafe" %in% methods && k >= 3L) {
failsafe <- safe(
metafor::fsn(x = model$yi, vi = model$vi, type = "Rosenberg"),
"Fail-safe N"
)
if (!is.null(failsafe)) add_test(
"Rosenberg fail-safe N",
as.character(failsafe$fsnum %||% ""),
.r4vn_meta_p(failsafe$pval %||% NA_real_, p_digits),
"Exploratory file-drawer analysis"
)
}
if ("selection" %in% methods && k >= 10L) {
selection <- safe(
metafor::selmodel(model, type = "stepfun", steps = c(.025, .05, .10, .50)),
"Selection model"
)
if (!is.null(selection)) {
b <- as.numeric(selection$beta[1L])
se <- as.numeric(selection$se[1L])
crit <- stats::qnorm(1 - (1 - ci) / 2)
adjusted[[length(adjusted) + 1L]] <- data.frame(
Method = "Step-function selection model",
`Adjusted effect (95% CI)` = .r4vn_meta_ci_text(
.r4vn_meta_back(b, config, ni),
.r4vn_meta_back(b - crit * se, config, ni),
.r4vn_meta_back(b + crit * se, config, ni), 3
),
`Imputed studies` = NA_integer_,
stringsAsFactors = FALSE, check.names = FALSE
)
}
}
if ("selection" %in% methods && k < 10L) {
notes <- c(notes, "Selection models were not fitted with fewer than 10 studies.")
}
table <- if (length(rows)) do.call(rbind, rows) else NULL
adjusted_table <- if (length(adjusted)) do.call(rbind, adjusted) else NULL
list(
status = if (length(notes)) paste(unique(notes), collapse = " ") else NULL,
table = table, adjusted = adjusted_table, methods = methods,
egger = egger, begg = begg, trimfill = trimfill,
failsafe = failsafe, selection = selection
)
}
.r4vn_meta_small_sample <- function(d, method, ci, config, ni, requested,
selected, digits = 2) {
methods <- unique(c("z", "knha", "adhoc", selected))
rows <- lapply(methods, function(test) {
fit <- tryCatch(
.r4vn_meta_model(d$yi, d$vi, method, ci = ci, test = test),
error = function(e) NULL
)
if (is.null(fit)) return(NULL)
pr <- .r4vn_meta_prediction(fit, config, ci, ni)
data.frame(
Inference = switch(test, z = "Normal approximation",
t = "t distribution",
knha = "Hartung-Knapp",
adhoc = "Modified Hartung-Knapp (adhoc)"),
Selected = ifelse(identical(test, selected), "Yes", "No"),
df = if (test == "z") NA_real_ else max(1, fit$k - fit$p),
`Effect (95% CI)` = .r4vn_meta_ci_text(
pr$estimate, pr$lower, pr$upper, digits
),
stringsAsFactors = FALSE, check.names = FALSE
)
})
rows <- Filter(Negate(is.null), rows)
comparison <- if (length(rows)) do.call(rbind, rows) else NULL
status <- if (nrow(d) <= 10L && identical(selected, "adhoc")) {
paste(
"Few-study inference is active. The modified Hartung-Knapp method",
"prevents the adjustment from producing smaller standard errors than",
"the conventional random-effects analysis."
)
} else if (nrow(d) <= 10L) {
paste0(
"Only ", nrow(d), " studies are available; inference method `",
selected, "` was selected. The sensitivity comparison is shown."
)
} else if (identical(requested, "auto")) {
"The normal approximation is used because more than 10 studies are available."
} else {
paste0("Inference method `", selected,
"` was explicitly selected; the sensitivity comparison is shown.")
}
list(
studies = nrow(d), requested = requested, selected = selected,
threshold = 10L,
status = status,
comparison = comparison
)
}
.r4vn_meta_leaveout <- function(model, config, digits, ni = NULL) {
if (model$k < 3L) return(NULL)
z <- tryCatch(metafor::leave1out(model), error = function(e) NULL)
if (is.null(z)) return(NULL)
zz <- as.data.frame(z)
if (!all(c("estimate", "ci.lb", "ci.ub") %in% names(zz))) return(NULL)
est <- .r4vn_meta_back(zz$estimate, config, ni)
lo <- .r4vn_meta_back(zz$ci.lb, config, ni)
hi <- .r4vn_meta_back(zz$ci.ub, config, ni)
slab <- rownames(zz)
if (is.null(slab) || !length(slab)) {
slab <- paste0("Study ", seq_len(nrow(zz)))
}
i2 <- if ("I2" %in% names(zz)) zz$I2 else rep(NA_real_, nrow(zz))
out <- data.frame(
`Study omitted` = slab,
`Effect (95% CI)` = mapply(
.r4vn_meta_ci_text, est, lo, hi,
MoreArgs = list(digits = digits), USE.NAMES = FALSE
),
`I-squared` = ifelse(
is.finite(i2), paste0(.r4vn_meta_fmt(i2, 1), "%"), ""
),
stringsAsFactors = FALSE, check.names = FALSE
)
list(raw = z, table = out, estimate = est, lower = lo, upper = hi)
}
.r4vn_meta_influence <- function(model, study) {
if (model$k < 3L) return(NULL)
inf <- tryCatch(stats::influence(model), error = function(e) NULL)
if (is.null(inf)) return(NULL)
d <- inf$inf
if (is.null(d)) return(list(raw = inf, table = NULL))
d <- as.data.frame(d)
pick <- function(candidates) {
hit <- intersect(candidates, names(d))
if (!length(hit)) rep(NA_real_, nrow(d)) else as.numeric(d[[hit[1L]]])
}
cook <- pick(c("cook.d", "cook.d."))
rstudent <- pick(c("rstudent", "rstudent."))
hat <- pick(c("hat", "hat."))
is_infl <- inf$is.infl
if (is.null(is_infl) || length(is_infl) != nrow(d)) {
is_infl <- rep(FALSE, nrow(d))
}
if (length(study) != nrow(d)) study <- rownames(d)
tab <- data.frame(
Study = study,
`Studentized residual` = vapply(
rstudent, .r4vn_meta_fmt, character(1), digits = 3
),
`Cook's distance` = vapply(
cook, .r4vn_meta_fmt, character(1), digits = 3
),
`Hat value` = vapply(
hat, .r4vn_meta_fmt, character(1), digits = 3
),
Influential = ifelse(is_infl, "Yes", "No"),
stringsAsFactors = FALSE, check.names = FALSE
)
list(raw = inf, table = tab, cook = cook, rstudent = rstudent,
hat = hat, influential = is_infl)
}
.r4vn_meta_cumulative <- function(d, order_name, method, test, ci, config,
digits, fixed = FALSE) {
if (is.null(order_name)) return(NULL)
order_label <- .r4vn_meta_variable_label(d[[order_name]], order_name)
ord <- order(d[[order_name]], na.last = NA)
if (length(ord) < 2L) return(NULL)
dd <- d[ord, , drop = FALSE]
rows <- vector("list", nrow(dd) - 1L)
for (i in 2:nrow(dd)) {
fit <- .r4vn_meta_model(
dd$yi[seq_len(i)], dd$vi[seq_len(i)], method, ci = ci, test = test,
fixed = fixed
)
pr <- .r4vn_meta_prediction(
fit, config, ci, dd$.n_for_back[seq_len(i)]
)
rows[[i - 1L]] <- data.frame(
`Up to` = as.character(dd[[order_name]][i]),
Studies = i,
`Effect (95% CI)` = .r4vn_meta_ci_text(
pr$estimate, pr$lower, pr$upper, digits
),
`I-squared` = paste0(.r4vn_meta_fmt(fit$I2, 1), "%"),
Estimate = pr$estimate,
Lower = pr$lower,
Upper = pr$upper,
stringsAsFactors = FALSE, check.names = FALSE
)
}
plot_data <- do.call(rbind, rows)
names(plot_data)[1L] <- paste0("Up to: ", order_label)
tab <- plot_data
tab$Estimate <- tab$Lower <- tab$Upper <- NULL
list(table = tab, plot_data = plot_data, order = order_name,
label = order_label)
}
.r4vn_meta_interpretation <- function(x) {
rows <- list()
add <- function(section, text) {
if (is.null(text) || !length(text) || is.na(text[1L]) ||
!nzchar(trimws(as.character(text)[1L]))) return(invisible(NULL))
rows[[length(rows) + 1L]] <<- data.frame(
Section = section,
Interpretation = as.character(text)[1L],
stringsAsFactors = FALSE, check.names = FALSE
)
invisible(NULL)
}
ptxt <- function(p) {
if (is.null(p) || !length(p) || !is.finite(p[1L])) return("not available")
paste0("p ", if (p[1L] < 10^(-x$p_digits)) {
paste0("< ", formatC(10^(-x$p_digits), format = "f",
digits = x$p_digits))
} else {
paste0("= ", formatC(p[1L], format = "f", digits = x$p_digits))
})
}
ci <- function(est, lo, hi, digits = x$digits) {
paste0(.r4vn_meta_fmt(est, digits), " (", round(x$ci * 100), "% CI ",
.r4vn_meta_fmt(lo, digits), " to ",
.r4vn_meta_fmt(hi, digits), ")")
}
null <- x$config$ref
crosses_null <- function(lo, hi) {
is.finite(lo) && is.finite(hi) && lo <= null && hi >= null
}
pr <- x$overall
h <- x$primary_model
model_name <- if (x$random) "random-effects" else "fixed-effect"
has_comparative_null <- !x$effect_type %in% c("PROP", "RATE")
direction <- if (!is.finite(pr$estimate) || !has_comparative_null) {
""
} else if (pr$estimate < null) {
paste0(" The pooled estimate was below the null value of ", null, ".")
} else if (pr$estimate > null) {
paste0(" The pooled estimate was above the null value of ", null, ".")
} else {
paste0(" The pooled estimate equaled the null value of ", null, ".")
}
null_text <- if (!has_comparative_null) {
""
} else if (crosses_null(pr$lower, pr$upper)) {
" The confidence interval included the null value."
} else {
" The confidence interval excluded the null value."
}
inference_name <- switch(
x$test_method, z = "normal-approximation", t = "t-distribution",
knha = "Hartung-Knapp", adhoc = "modified Hartung-Knapp",
x$test_method
)
add(
"Overall effect",
paste0(
h$k, " studies were included. The ", model_name,
" model estimated a pooled ", x$config$label, " of ",
ci(pr$estimate, pr$lower, pr$upper), "; ", ptxt(pr$p), ". ",
if (is.finite(pr$p) && pr$p < (1 - x$ci)) {
"The pooled estimate was statistically significant."
} else if (is.finite(pr$p)) {
"The pooled estimate was not statistically significant."
} else {
"Statistical significance of the pooled estimate was not available."
},
direction, null_text, " Inference used the ", inference_name,
" method."
)
)
i2 <- as.numeric(h$I2)
i2_description <- if (!is.finite(i2)) {
"not estimable"
} else if (i2 < 25) {
"low"
} else if (i2 < 50) {
"moderate"
} else if (i2 < 75) {
"substantial"
} else {
"considerable"
}
add(
"Heterogeneity",
paste0(
"Between-study heterogeneity was ", i2_description,
" on the I-squared scale (I-squared = ",
.r4vn_meta_fmt(i2, 1), "%; tau-squared = ",
.r4vn_meta_fmt(h$tau2, 4), "). Cochran's Q was ",
.r4vn_meta_fmt(h$QE, 2), " with ", max(0, h$k - h$p),
" degrees of freedom (", ptxt(h$QEp), "). ",
if (is.finite(h$QEp) && h$QEp < (1 - x$ci)) {
"Heterogeneity was statistically significant."
} else if (is.finite(h$QEp)) {
"There was no statistically significant evidence of heterogeneity."
} else {
"Statistical significance of heterogeneity was not available."
}
)
)
if (isTRUE(x$prediction) && isTRUE(x$random) &&
is.finite(pr$pi_lower) && is.finite(pr$pi_upper)) {
pi_text <- if (!has_comparative_null) {
NULL
} else if (crosses_null(pr$pi_lower, pr$pi_upper)) {
"included"
} else {
"excluded"
}
add(
"Prediction interval",
paste0(
"The ", round(x$ci * 100), "% prediction interval was ",
.r4vn_meta_fmt(pr$pi_lower, x$digits), " to ",
.r4vn_meta_fmt(pr$pi_upper, x$digits),
if (is.null(pi_text)) "." else paste0(" and ", pi_text, " the null value."),
" This interval describes the range expected for the ",
"underlying effect in a comparable new study, conditional on the fitted model."
)
)
}
if (length(x$subgroups)) {
for (nm in names(x$subgroups)) {
z <- x$subgroups[[nm]]
fit <- z$test_model
if (is.null(fit)) next
evidence <- if (is.finite(fit$QMp) && fit$QMp < 0.05) {
"There was statistical evidence that pooled effects differed between subgroups."
} else {
"There was no statistical evidence that pooled effects differed between subgroups."
}
add(
paste0("Subgroup analysis: ", z$label %||% nm),
paste0(
evidence, " The test for subgroup differences gave Q = ",
.r4vn_meta_fmt(fit$QM, 2), " with ", fit$m,
" degrees of freedom (", ptxt(fit$QMp), "). This test should be ",
"interpreted cautiously when subgroups contain few studies."
)
)
}
}
if (!is.null(x$regression$model)) {
fit <- x$regression$model
sm <- summary(fit)
beta_names <- rownames(sm$beta)
if (is.null(beta_names)) beta_names <- paste0("coefficient ", seq_along(sm$pval))
display_beta_names <- .r4vn_meta_term_labels(
beta_names, x$analysis_data, x$reg
)
moderator_idx <- !grepl("^(intrcpt|intercept|\\(intercept\\))$",
beta_names, ignore.case = TRUE)
sig <- display_beta_names[
moderator_idx & is.finite(sm$pval) & sm$pval < 0.05
]
sig_text <- if (length(sig)) {
paste0(" Statistically associated coefficient(s): ",
paste(sig, collapse = ", "), ".")
} else {
" No individual moderator coefficient met the 0.05 significance threshold."
}
add(
"Meta-regression",
paste0(
"The joint moderator test gave QM = ", .r4vn_meta_fmt(fit$QM, 2),
" with ", fit$m, " degrees of freedom (", ptxt(fit$QMp), ").",
sig_text, " Residual heterogeneity was I-squared = ",
.r4vn_meta_fmt(fit$I2, 1), "% and tau-squared = ",
.r4vn_meta_fmt(fit$tau2, 4),
". Meta-regression findings are observational across studies and should not be interpreted causally."
)
)
}
if (!is.null(x$small_sample$status)) {
add("Few-study inference", x$small_sample$status)
}
if (!is.null(x$bias)) {
bias_bits <- character()
if (!is.null(x$bias$egger)) {
bias_bits <- c(
bias_bits,
paste0(
"Egger's regression test ",
if (x$bias$egger$pval < 0.05) "detected" else "did not detect",
" funnel asymmetry (", ptxt(x$bias$egger$pval), ")"
)
)
}
if (!is.null(x$bias$begg)) {
bias_bits <- c(
bias_bits,
paste0(
"the rank-correlation test ",
if (x$bias$begg$pval < 0.05) "detected" else "did not detect",
" asymmetry (", ptxt(x$bias$begg$pval), ")"
)
)
}
if (!is.null(x$bias$trimfill)) {
adj <- .r4vn_meta_prediction(
x$bias$trimfill, x$config, x$ci, x$analysis_data$.n_for_back
)
bias_bits <- c(
bias_bits,
paste0(
"trim-and-fill imputed ", x$bias$trimfill$k0 %||% 0L,
" study/studies and gave an adjusted pooled estimate of ",
ci(adj$estimate, adj$lower, adj$upper)
)
)
}
if (length(bias_bits)) {
add(
"Small-study effects / publication bias",
paste0(
paste(bias_bits, collapse = "; "),
". Funnel asymmetry and trim-and-fill are sensitivity diagnostics and do not by themselves prove or exclude publication bias."
)
)
} else if (!is.null(x$bias$status)) {
add("Small-study effects / publication bias", x$bias$status)
}
}
if (!is.null(x$influence$table)) {
influential <- x$influence$table$Study[
x$influence$table$Influential %in% "Yes"
]
add(
"Influence diagnostics",
if (length(influential)) {
paste0(
"The following studies were flagged as influential: ",
paste(influential, collapse = ", "),
". Their data and analytic assumptions should be reviewed."
)
} else {
"No study was flagged as influential by the fitted influence diagnostics."
}
)
}
if (!is.null(x$leaveout$estimate) && length(x$leaveout$estimate)) {
rng <- range(x$leaveout$estimate, finite = TRUE)
add(
"Leave-one-out analysis",
paste0(
"Across leave-one-out analyses, pooled estimates ranged from ",
.r4vn_meta_fmt(rng[1L], x$digits), " to ",
.r4vn_meta_fmt(rng[2L], x$digits),
". Compare this range with the main pooled estimate to assess robustness to individual studies."
)
)
}
if (!is.null(x$cumulative$plot_data) && nrow(x$cumulative$plot_data)) {
z <- x$cumulative$plot_data
add(
"Cumulative meta-analysis",
paste0(
"The cumulative pooled estimate changed from ",
ci(z$Estimate[1L], z$Lower[1L], z$Upper[1L]), " at ",
x$cumulative$label %||% x$cumulative$order, " = ", z[[1L]][1L],
" to ", ci(z$Estimate[nrow(z)], z$Lower[nrow(z)], z$Upper[nrow(z)]),
" at ", x$cumulative$label %||% x$cumulative$order, " = ",
z[[1L]][nrow(z)],
". The cumulative table and plot show how evidence evolved as studies accumulated."
)
)
}
if (!length(rows)) {
return(data.frame(Section = character(), Interpretation = character(),
stringsAsFactors = FALSE))
}
out <- do.call(rbind, rows)
rownames(out) <- NULL
out
}
.r4vn_meta_results_text <- function(x) {
z <- .r4vn_meta_interpretation(x)
paste(z$Interpretation, collapse = " ")
}
.r4vn_meta_html_table <- function(data, title = NULL) {
if (is.null(data) || !is.data.frame(data)) return("")
th <- paste0("<th>", .r4vn_meta_escape(names(data)), "</th>",
collapse = "")
body <- if (!nrow(data)) "" else paste(
vapply(seq_len(nrow(data)), function(i) {
vals <- vapply(data[i, , drop = FALSE], function(z) {
if (is.na(z)) "" else as.character(z)
}, character(1))
paste0(
"<tr>",
paste0("<td>", .r4vn_meta_escape(vals), "</td>", collapse = ""),
"</tr>"
)
}, character(1)),
collapse = ""
)
heading <- if (is.null(title) || !nzchar(title)) "" else
paste0("<h2>", .r4vn_meta_escape(title), "</h2>")
paste0(
heading,
"<div class=\"table-wrap\"><table><thead><tr>", th,
"</tr></thead><tbody>", body, "</tbody></table></div>"
)
}
.r4vn_meta_report_html <- function(x) {
sections <- character()
sections <- c(sections, .r4vn_meta_html_table(x$overview, "Analysis overview"))
sections <- c(sections, .r4vn_meta_html_table(x$table, "Main analysis"))
sections <- c(
sections,
.r4vn_meta_html_table(x$significance_tests, "Statistical significance tests")
)
sections <- c(
sections,
.r4vn_meta_html_table(x$heterogeneity, "Model and heterogeneity")
)
if (length(x$subgroup_tables)) {
for (nm in names(x$subgroup_tables)) {
subgroup_label <- x$subgroups[[nm]]$label %||% nm
sections <- c(
sections,
.r4vn_meta_html_table(
x$subgroup_tables[[nm]],
paste0("Subgroup analysis: ", subgroup_label)
)
)
}
sections <- c(
sections,
.r4vn_meta_html_table(
x$tables$Subgroup_test, "Test for subgroup differences"
)
)
}
if (!is.null(x$regression$table)) {
sections <- c(
sections,
.r4vn_meta_html_table(x$regression$table, "Meta-regression"),
.r4vn_meta_html_table(x$regression$statistics,
"Meta-regression statistics")
)
if (!is.null(x$regression$univariable)) {
sections <- c(
sections,
.r4vn_meta_html_table(x$regression$univariable,
"Univariable moderator screening")
)
}
}
if (!is.null(x$bias$table)) {
sections <- c(
sections,
.r4vn_meta_html_table(
x$bias$table, "Small-study effects / publication bias"
)
)
}
if (!is.null(x$bias$status)) {
sections <- c(
sections,
paste0("<p class=\"note\">", .r4vn_meta_escape(x$bias$status), "</p>")
)
}
if (!is.null(x$bias$adjusted)) {
sections <- c(
sections,
.r4vn_meta_html_table(x$bias$adjusted,
"Publication-bias sensitivity estimates")
)
}
if (!is.null(x$small_sample$comparison)) {
sections <- c(
sections,
.r4vn_meta_html_table(x$small_sample$comparison,
"Few-study inference sensitivity"),
paste0("<p class=\"note\">",
.r4vn_meta_escape(x$small_sample$status), "</p>")
)
}
if (!is.null(x$leaveout$table)) {
sections <- c(
sections,
.r4vn_meta_html_table(
x$leaveout$table, "Leave-one-out sensitivity analysis"
)
)
}
if (!is.null(x$influence$table)) {
sections <- c(
sections,
.r4vn_meta_html_table(x$influence$table, "Influence diagnostics")
)
}
if (!is.null(x$cumulative$table)) {
sections <- c(
sections,
.r4vn_meta_html_table(
x$cumulative$table, "Cumulative meta-analysis"
)
)
}
if (isTRUE(x$report)) {
narrative <- if (!is.null(x$interpretation) && nrow(x$interpretation)) {
paste(
vapply(seq_len(nrow(x$interpretation)), function(i) {
paste0(
"<h3>", .r4vn_meta_escape(x$interpretation$Section[i]), "</h3>",
"<p>", .r4vn_meta_escape(x$interpretation$Interpretation[i]), "</p>"
)
}, character(1)),
collapse = "\n"
)
} else {
paste0("<p>", .r4vn_meta_escape(x$results_text), "</p>")
}
sections <- c(sections, paste0("<h2>Interpretation</h2>", narrative))
}
if (length(x$plots)) {
for (nm in names(x$plots)) {
sections <- c(
sections,
paste0(
"<h2>", .r4vn_meta_escape(x$plot_titles[[nm]]), "</h2>",
"<div class=\"figure\"><img src=\"",
.r4vn_meta_escape(basename(x$plots[[nm]])),
"\" alt=\"", .r4vn_meta_escape(x$plot_titles[[nm]]),
"\"></div>"
)
)
}
}
css <- paste0(
"body{font-family:'Times New Roman',Times,serif;color:#111;background:#fff;",
"margin:24px;line-height:1.35}",
".report{max-width:1180px;margin:auto}",
"h1{font-size:22px;margin:0 0 16px}",
"h2{font-size:18px;margin:24px 0 8px}",
"h3{font-size:15px;margin:15px 0 5px}",
".table-wrap{overflow-x:auto;margin-bottom:14px}",
"table{border-collapse:collapse;width:100%;border-top:2px solid #111;",
"border-bottom:2px solid #111}",
"th{border-bottom:1.5px solid #111;padding:6px 8px;text-align:center;",
"font-weight:700;white-space:nowrap}",
"td{padding:5px 8px;text-align:center;vertical-align:top}",
"th:first-child,td:first-child{text-align:left}",
".note{font-size:13px}",
".figure{text-align:center;margin:8px 0 22px}",
".figure img{max-width:100%;height:auto}",
".foot{font-size:12px;margin-top:18px;color:#444}"
)
paste0(
"<!doctype html><html><head><meta charset=\"utf-8\">",
"<meta name=\"viewport\" content=\"width=device-width,initial-scale=1\">",
"<title>", .r4vn_meta_escape(x$title), "</title><style>", css,
"</style></head><body><main class=\"report\"><h1>",
.r4vn_meta_escape(x$title), "</h1>",
paste(sections, collapse = "\n"),
"<div class=\"foot\">R4VN meta-analysis. Random-effects estimator: ",
.r4vn_meta_escape(x$method), "; inference: ",
.r4vn_meta_escape(x$test_method), ".",
"</div></main></body></html>"
)
}
.r4vn_meta_open <- function(path) {
viewer <- getOption("viewer")
p <- normalizePath(path, winslash = "/", mustWork = TRUE)
if (is.function(viewer)) viewer(p) else utils::browseURL(p)
invisible(path)
}
#' Publication-Ready Meta-Analysis in One Command
#'
#' Fits fixed- and/or random-effects meta-analysis and returns a complete,
#' publication-ready report. `tabmeta()` accepts a binary 2-by-2 table
#' (`a`, `b`, `c`, `d`), event/total data, continuous summaries, rates,
#' correlations, or study-level effect estimates with standard errors or
#' confidence intervals. The default `profile="auto"` adds prediction,
#' few-study inference, subgroup/moderator results when requested,
#' small-study-effect diagnostics, sensitivity analyses, and figures.
#'
#' @param data Optional data frame. When omitted, active R4VN data are used.
#' @param study Study label variable.
#' @param effect Generic study-level effect estimate on the natural scale.
#' @param se Standard error on the analysis scale. For ratio measures this is
#' the standard error of the log effect.
#' @param lower,upper Lower and upper confidence limits for `effect`.
#' @param a,b,c,d Binary 2-by-2 cells: events and non-events in group 1
#' (`a`, `b`) and group 0 (`c`, `d`). Supplying these four arguments is
#' equivalent to `event1=a`, `n1=a+b`, `event0=c`, `n0=c+d`.
#' @param event1,n1,event0,n0 Events and total sample sizes in groups 1 and 0.
#' @param mean1,sd1,mean0,sd0 Group means and standard deviations.
#' @param event,n,time Single-group events, sample size, and person-time.
#' @param time1,time0 Person-time in groups 1 and 0 for incidence-rate ratios.
#' @param or,rr,rd,hr,irr,md,smd,prop,rate,cor Logical effect selectors.
#' With raw binary input, OR is inferred if none is selected; with continuous
#' summaries, MD is inferred. Set a selector to request another measure.
#' @param fixed Fit a fixed-effect model in addition to, or instead of, the
#' random-effects model.
#' @param random Fit a random-effects model; default `TRUE`.
#' @param method Random-effects tau-squared estimator; default `"REML"`.
#' @param hk Backward-compatible logical shortcut. `TRUE` is
#' `small="knha"`; `FALSE` is `small="z"`. Prefer `small` in new code.
#' @param small Random-effects inference: `"auto"`, `"adhoc"`, `"knha"`,
#' `"t"`, or `"z"`. `"auto"` uses modified Hartung-Knapp (`"adhoc"`)
#' when there are at most 10 studies and the normal approximation otherwise.
#' @param prediction Add a prediction interval. `NULL` uses the profile default.
#' @param by,subgroup One or more subgroup variables. Use one unquoted variable,
#' a character vector, or `vars(region, design)`. `subgroup` is a readable
#' alias for `by`; both may be combined. Factor levels and value labels from
#' labelled data are used in tables, interpretation, and subgroup figures.
#' For every estimable subgroup, the result includes its pooled effect and
#' confidence interval, effect test/df/p-value, prediction interval,
#' Cochran's Q/df/p-value, I-squared and tau-squared with confidence
#' intervals when estimable. The publication table is transposed: statistics
#' are rows and subgroup labels are columns. Multiple subgroup variables
#' produce one transposed table per variable.
#' @param reg,moderator One or more meta-regression moderators created with
#' `vars()` or supplied as character names. `moderator` is an alias for
#' `reg`. The output includes the multivariable model and a univariable
#' moderator screen.
#' @param bias Assess small-study effects/publication bias. `NULL` uses the
#' profile default.
#' @param bias_methods One or more of `"egger"`, `"begg"`, `"trimfill"`,
#' `"failsafe"`, or `"selection"`. `"auto"` runs Egger, Begg, and
#' trim-and-fill; the `"full"` profile also requests fail-safe N and a
#' selection model. Methods are used as sensitivity diagnostics, not as proof
#' that publication bias is or is not present.
#' @param leaveout Perform leave-one-out sensitivity analysis. `NULL` uses the
#' profile default.
#' @param influence Perform influence diagnostics. `NULL` uses the profile
#' default.
#' @param cumulative Optional variable defining the ordering for cumulative
#' meta-analysis, usually publication year.
#' @param transform Single-proportion transformation: `"logit"`, `"arcsine"`,
#' `"ft"`, or `"none"`.
#' @param cc Continuity correction for zero events; default `0.5`.
#' @param zero Handling of double-zero binary studies: `"keep"` or `"exclude"`.
#' @param ci Confidence level as a proportion.
#' @param digit Number of decimals for effect estimates and confidence limits,
#' including subgroup and meta-regression estimates. Default `2`.
#' @param p_digit Number of decimals for p-values.
#' @param profile Analysis profile. `"auto"` (default) creates a complete
#' appropriate analysis; `"brief"` creates the main model and forest plot;
#' `"full"` also requests extended bias diagnostics and radial/L'Abbe plots;
#' `"custom"` turns optional modules off unless explicitly requested.
#' @param full Backward-compatible profile shortcut: `TRUE` selects `"full"`
#' and `FALSE` selects `"custom"`.
#' @param plot Create publication-ready plots. `NULL` uses the profile default.
#' With subgroup analysis, the HTML Viewer includes the overall forest plot
#' and one clearly titled forest plot for every estimable subgroup level.
#' @param plot_display Plot type(s) also drawn in the interactive R/RStudio Plot
#' pane when `plot=TRUE`. Default `"forest"`. Use a character vector for Plot
#' history, `"all"` for every available figure, or `"none"` to suppress Plot
#' pane drawing while retaining figures in the HTML Viewer. When subgroup
#' analysis is present, selecting `"forest"` adds the overall forest plot and
#' every subgroup forest plot to Plot history; use Previous/Next to review.
#' @param plot_args Named list of plot options. Supply common options directly,
#' or nested lists such as `list(all=list(color="navy"),
#' forest=list(xlim=c(0.2, 2)))`.
#' @param report Add manuscript-style results text.
#' @param interpretation Add a detailed, sectioned interpretation covering the
#' pooled effect, heterogeneity, prediction interval, subgroup differences,
#' moderators, small-study effects, few-study inference, influence,
#' leave-one-out robustness, and cumulative evidence whenever available.
#' Default is `FALSE`.
#' @param title Report title.
#' @param show Open the complete HTML report; default `TRUE`.
#' @param export Optional direct export format(s): `"html"`, `"docx"`,
#' `"xlsx"`, `"pdf"`, or `"png"`.
#' @param file Export file or base path. The extension may determine the format.
#' @param open Open the exported file.
#' @param strict Stop when an optional diagnostic or figure cannot be created.
#' The default `FALSE` retains the main analysis and issues a warning/note.
#'
#' @return Invisibly returns an object of class `r4vn_meta` and `r4vn_tab`.
#' Stable publication components are available in `$estimates`, `$tests`,
#' `$diagnostics`, `$models`, `$tables`, and `$metadata`. The table
#' `$tables$Statistical_tests` gives the pooled-effect and Cochran's Q tests.
#' Its significance and conclusion columns are added only when
#' `interpretation=TRUE`. `$tables$Subgroup` is the transposed publication
#' table for the first subgroup variable, `$subgroup_tables` contains every
#' transposed subgroup table, and `$subgroup_long` retains tidy long output.
#' @family R4VN tables
#' @seealso `tabexport`, `vars`
#' @export
#'
#' @examples
#' if (requireNamespace("metafor", quietly = TRUE)) {
#' dat <- read.csv(
#' system.file("extdata", "meta_example.csv", package = "R4VN")
#' )
#'
#' # 1. Binary outcome from a, b, c, d. OR is inferred automatically.
#' dat$non_event_treat <- dat$n_treat - dat$event_treat
#' dat$non_event_control <- dat$n_control - dat$event_control
#' m_abcd <- tabmeta(
#' data = dat, study = study,
#' a = event_treat, b = non_event_treat,
#' c = event_control, d = non_event_control,
#' profile = "custom", plot = FALSE, show = FALSE
#' )
#' m_abcd$estimates$overall
#' m_abcd$tables$Binary_2x2
#' m_abcd$tables$Statistical_tests
#' m_abcd$tests$overall
#' m_abcd$tests$heterogeneity
#' plot(
#' m_abcd, type = "forest", show_abcd = TRUE,
#' abcd_titles = c("Events T", "No event T", "Events C", "No event C")
#' )
#'
#' # 2. Equivalent event/total syntax; request RR or RD with rr/rd=TRUE.
#' m_or <- tabmeta(
#' data = dat, study = study,
#' event1 = event_treat, n1 = n_treat,
#' event0 = event_control, n0 = n_control,
#' or = TRUE, profile = "custom", plot = FALSE, show = FALSE
#' )
#'
#' # 3. Complete automatic analysis: report, prediction, diagnostics, plots.
#' \donttest{
#' m_all <- tabmeta(
#' data = dat, study = study,
#' a = event_treat, b = non_event_treat,
#' c = event_control, d = non_event_control,
#' show = FALSE
#' )
#' }
#'
#' # 4. Subgroup analysis: one or several subgroup variables.
#' m_sub <- tabmeta(
#' data = dat, study = study,
#' effect = OR, lower = LCI, upper = UCI, or = TRUE,
#' subgroup = region,
#' profile = "custom", plot = FALSE, show = FALSE
#' )
#' m_sub$tables$Subgroup
#' m_sub$subgroup_long
#' m_sub$tables$Subgroup_test
#' # Each subgroup row includes its effect test and heterogeneity test.
#' names(m_sub$tables$Subgroup)
#' # If the subgroup variable carries value labels, R4VN prints those labels
#' # instead of numeric codes. With plot=TRUE, both Viewer and Plot history
#' # contain the overall plot and one clearly titled plot per subgroup level.
#' dat$risk_group <- rep(c(1, 2), length.out = nrow(dat))
#' attr(dat$risk_group, "label") <- "Baseline risk"
#' attr(dat$risk_group, "labels") <- c("Lower risk" = 1, "Higher risk" = 2)
#' \donttest{
#' m_labelled <- tabmeta(
#' data = dat, study = study,
#' effect = OR, lower = LCI, upper = UCI, or = TRUE,
#' subgroup = risk_group, profile = "custom",
#' plot = TRUE, plot_display = "forest", show = FALSE
#' )
#' m_labelled$plot_titles
#' }
#'
#' # Multiple subgroup analyses can be requested together.
#' dat$period <- ifelse(dat$year < median(dat$year), "Earlier", "Later")
#' m_sub2 <- tabmeta(
#' data = dat, study = study,
#' effect = OR, lower = LCI, upper = UCI, or = TRUE,
#' subgroup = vars(region, period),
#' profile = "custom", plot = FALSE, show = FALSE
#' )
#'
#' # 5. Meta-regression with numeric and categorical moderators.
#' attr(dat$year, "label") <- "Publication year"
#' attr(dat$region, "label") <- "Geographic region"
#' m_reg <- tabmeta(
#' data = dat, study = study,
#' effect = OR, lower = LCI, upper = UCI, or = TRUE,
#' moderator = vars(c.year, region),
#' profile = "custom", plot = FALSE, show = FALSE
#' )
#' m_reg$tables$Meta_regression
#' m_reg$tables$Moderator_univariable
#' # Intercept is written in full; moderator labels are used when available.
#' # Set digit=4, for example, when four decimal places are required.
#'
#' # 6. Publication-bias and small-study-effect sensitivity analyses.
#' m_bias <- tabmeta(
#' data = dat, study = study,
#' effect = OR, lower = LCI, upper = UCI, or = TRUE,
#' bias = TRUE,
#' bias_methods = c("egger", "begg", "trimfill", "failsafe"),
#' profile = "custom", plot = FALSE, show = FALSE
#' )
#' m_bias$tables$Publication_bias
#' m_bias$tables$Publication_bias_adjusted
#'
#' # 7. Few-study inference. auto uses modified Hartung-Knapp at <=10 studies.
#' m_few <- tabmeta(
#' data = dat[1:8, ], study = study,
#' effect = OR, lower = LCI, upper = UCI, or = TRUE,
#' small = "auto", profile = "custom", plot = FALSE, show = FALSE
#' )
#' m_few$tables$Small_sample_inference
#'
#' # 8. Generic hazard ratios with confidence intervals.
#' m_hr <- tabmeta(
#' data = dat, study = study,
#' effect = OR, lower = LCI, upper = UCI, hr = TRUE,
#' profile = "custom", plot = FALSE, show = FALSE
#' )
#'
#' # 9. Continuous MD/SMD, single proportion/rate, correlation, and IRR.
#' cont <- data.frame(
#' study = paste0("C", 1:5),
#' m1 = c(12, 14, 13, 16, 15), s1 = c(3, 4, 3, 5, 4), n1 = rep(60, 5),
#' m0 = c(15, 15, 16, 18, 16), s0 = c(4, 4, 5, 5, 4), n0 = rep(60, 5)
#' )
#' m_md <- tabmeta(
#' cont, study, mean1 = m1, sd1 = s1, n1 = n1,
#' mean0 = m0, sd0 = s0, n0 = n0, md = TRUE,
#' profile = "custom", plot = FALSE, show = FALSE
#' )
#' m_smd <- tabmeta(
#' cont, study, mean1 = m1, sd1 = s1, n1 = n1,
#' mean0 = m0, sd0 = s0, n0 = n0, smd = TRUE,
#' profile = "custom", plot = FALSE, show = FALSE
#' )
#'
#' one <- data.frame(
#' study = paste0("P", 1:5), events = c(8, 12, 15, 10, 14),
#' total = c(100, 110, 120, 90, 105), person_time = c(80, 90, 95, 75, 88),
#' correlation = c(.20, .28, .15, .31, .24)
#' )
#' m_prop <- tabmeta(
#' one, study, event = events, n = total, prop = TRUE,
#' profile = "custom", plot = FALSE, show = FALSE
#' )
#' m_rate <- tabmeta(
#' one, study, event = events, time = person_time, rate = TRUE,
#' profile = "custom", plot = FALSE, show = FALSE
#' )
#' m_cor <- tabmeta(
#' one, study, effect = correlation, n = total, cor = TRUE,
#' profile = "custom", plot = FALSE, show = FALSE
#' )
#' m_irr <- tabmeta(
#' dat[1:5, ], study,
#' event1 = event_treat, time1 = n_treat,
#' event0 = event_control, time0 = n_control, irr = TRUE,
#' profile = "custom", plot = FALSE, show = FALSE
#' )
#'
#' # 10. Cumulative meta-analysis ordered by publication year.
#' m_cum <- tabmeta(
#' data = dat, study = study,
#' effect = OR, lower = LCI, upper = UCI, or = TRUE,
#' cumulative = year,
#' profile = "custom", plot = FALSE, show = FALSE
#' )
#' m_cum$tables$Cumulative
#'
#' # 11. Draw or save individual publication figures.
#' if (interactive()) {
#' # plot=TRUE keeps all figures in the HTML Viewer and also draws the
#' # selected figures in the R/RStudio Plot pane and Plot history.
#' m_publication <- tabmeta(
#' data = dat, study = study,
#' a = event_treat, b = non_event_treat,
#' c = event_control, d = non_event_control,
#' interpretation = TRUE,
#' plot = TRUE,
#' plot_display = c("forest", "funnel", "trimfill"),
#' plot_args = list(
#' all = list(
#' font_family = "Arial", background = "white",
#' title_color = "#17365D", title_size = 1.1,
#' text_size = 0.86, axis_size = 0.92
#' ),
#' forest = list(
#' subtitle = "Random-effects model with 95% confidence intervals",
#' caption = "Square size reflects study weight; diamond is pooled effect.",
#' margins = c(5.5, 4.2, 5.0, 2.0),
#' point_color = "#1F4E79", ci_color = "#5B9BD5",
#' summary_color = "#C00000", summary_border = "#7F0000",
#' point_shape = 15, row_shade = "zebra",
#' shade_color = "#F5F7FA", show_weights = TRUE,
#' weight_title = "Weight", estimate_title = "OR (95% CI)",
#' show_prediction = FALSE,
#' ref_color = "#666666", ref_type = 2,
#' xlim = c(0.2, 2.0), ticks = c(0.25, 0.5, 1, 1.5, 2)
#' ),
#' funnel = list(
#' point_shape = 21, point_color = "#1F4E79",
#' point_bg = "#D9EAF7", point_size = 1.1,
#' contour_levels = c(90, 95, 99),
#' contour_colors = c("#FFF2CC", "#FCE4D6", "#E2F0D9"),
#' funnel_label = "out", funnel_legend = "topright"
#' )
#' )
#' )
#' m_publication$tables$Interpretation
#'
#' # Any figure can be redrawn or saved independently.
#' plot(m_publication, type = "forest")
#' plot(m_publication, type = "funnel", contour = TRUE)
#' plot(m_publication, type = "trimfill", contour = TRUE)
#' plot(m_reg, type = "bubble", moderator = "year")
#' plot(
#' m_publication, type = "forest", file = "forest_publication.tiff",
#' width = 2400, height = 1800, res = 300,
#' font_family = "Arial", point_color = "#1F4E79",
#' ci_color = "#5B9BD5", summary_color = "#C00000",
#' show_weights = TRUE, show_prediction = FALSE,
#' xlim = c(0.5, 1.5), ticks = c(0.5, 0.75, 1, 1.25, 1.5)
#' )
#' # Values outside xlim remain exact in the Estimate (95% CI) column;
#' # the graphical confidence interval is clipped with an arrow.
#' m_publication$overall[c("pi_lower", "pi_upper")]
#' tabmeta(
#' data = dat, study = study,
#' a = event_treat, b = non_event_treat,
#' c = event_control, d = non_event_control,
#' export = c("docx", "xlsx"), file = "meta_report",
#' show = FALSE
#' )
#' }
#' }
tabmeta <- function(
data = NULL,
study,
effect = NULL, se = NULL, lower = NULL, upper = NULL,
a = NULL, b = NULL, c = NULL, d = NULL,
event1 = NULL, n1 = NULL, event0 = NULL, n0 = NULL,
mean1 = NULL, sd1 = NULL, mean0 = NULL, sd0 = NULL,
event = NULL, n = NULL, time = NULL, time1 = NULL, time0 = NULL,
or = FALSE, rr = FALSE, rd = FALSE, hr = FALSE, irr = FALSE,
md = FALSE, smd = FALSE, prop = FALSE, rate = FALSE, cor = FALSE,
fixed = FALSE, random = TRUE, method = "REML", hk = NULL,
small = c("auto", "adhoc", "knha", "t", "z"),
prediction = NULL, by = NULL, subgroup = NULL,
reg = NULL, moderator = NULL,
bias = NULL, bias_methods = "auto",
leaveout = NULL, influence = NULL, cumulative = NULL,
transform = NULL, cc = 0.5, zero = c("keep", "exclude"),
ci = 0.95, digit = 2, p_digit = 3,
profile = c("auto", "brief", "full", "custom"),
full = NULL, plot = NULL, plot_display = "forest", plot_args = list(),
report = FALSE, interpretation = FALSE,
title = NULL, show = TRUE,
export = NULL, file = NULL, open = FALSE, strict = FALSE) {
# `c` is a public 2-by-2 cell argument. Capture its expression before
# restoring base::c locally; otherwise the argument promise can mask c()
# and be forced while ordinary vectors are being constructed.
c_expr <- substitute(c)
c <- base::c
.r4vn_meta_need()
data <- .r4vn_resolve_analysis_data(data)
env <- parent.frame()
settings <- .r4vn_meta_profile(
profile = match.arg(profile), full = full,
prediction = prediction, bias = bias, leaveout = leaveout,
influence = influence, plot = plot
)
profile <- settings$profile
prediction <- settings$prediction
bias <- settings$bias
leaveout <- settings$leaveout
influence <- settings$influence
plot <- settings$plot
full <- settings$full
flag_names <- c(
"or", "rr", "rd", "hr", "irr", "md", "smd", "prop", "rate", "cor",
"fixed", "random", "prediction", "bias", "leaveout",
"influence", "full", "plot", "report", "interpretation", "show",
"open", "strict"
)
for (z in flag_names) .r4vn_meta_flag(get(z), z)
if (!is.null(hk)) .r4vn_meta_flag(hk, "hk")
if (!is.list(plot_args) ||
(length(plot_args) && is.null(names(plot_args)))) {
stop("`plot_args` must be a named list.", call. = FALSE)
}
plot_types <- c(
"forest", "funnel", "trimfill", "leaveout", "influence", "baujat",
"radial", "labbe", "cumulative", "bubble"
)
if (is.null(plot_display) || identical(plot_display, FALSE)) {
plot_display <- character()
} else {
plot_display <- unique(tolower(as.character(plot_display)))
plot_display <- plot_display[!is.na(plot_display) & nzchar(plot_display)]
if ("none" %in% plot_display) plot_display <- character()
bad_display <- setdiff(plot_display, c(plot_types, "all"))
if (length(bad_display)) {
stop("Unsupported `plot_display`: ", paste(bad_display, collapse = ", "),
".", call. = FALSE)
}
}
if (!fixed && !random) {
stop("At least one of `fixed` or `random` must be TRUE.", call. = FALSE)
}
.r4vn_meta_num1(ci, "ci", 0, 1, FALSE, FALSE)
.r4vn_meta_num1(cc, "cc", 0, Inf)
digit <- as.integer(.r4vn_meta_num1(digit, "digit", 0, 10))
p_digit <- as.integer(.r4vn_meta_num1(p_digit, "p_digit", 0, 10))
method_input <- toupper(as.character(method)[1L])
method_map <- c(
REML = "REML", ML = "ML", DL = "DL", HE = "HE", HS = "HS",
HSK = "HSk", SJ = "SJ", EB = "EB", PM = "PM", PMM = "PMM",
GENQ = "GENQ"
)
if (!method_input %in% names(method_map)) {
stop("Unsupported random-effects method: ", method_input, ".",
call. = FALSE)
}
method <- unname(method_map[[method_input]])
zero <- match.arg(zero)
small <- match.arg(tolower(as.character(small)[1L]),
c("auto", "adhoc", "knha", "t", "z"))
if (!is.null(hk)) small <- if (isTRUE(hk)) "knha" else "z"
bias_methods <- .r4vn_meta_bias_methods(bias_methods, profile)
report_text <- isTRUE(report) || isTRUE(interpretation)
expr <- list(
study = substitute(study),
effect = substitute(effect), se = substitute(se),
lower = substitute(lower), upper = substitute(upper),
a = substitute(a), b = substitute(b), c = c_expr, d = substitute(d),
event1 = substitute(event1), n1 = substitute(n1),
event0 = substitute(event0), n0 = substitute(n0),
mean1 = substitute(mean1), sd1 = substitute(sd1),
mean0 = substitute(mean0), sd0 = substitute(sd0),
event = substitute(event), n = substitute(n),
time = substitute(time), time1 = substitute(time1), time0 = substitute(time0),
by = substitute(by), subgroup = substitute(subgroup),
cumulative = substitute(cumulative)
)
study_value <- .r4vn_meta_eval(expr$study, data, env, "study", TRUE)
if (length(study_value) != nrow(data)) {
stop("`study` must have one value per row of `data`.", call. = FALSE)
}
study_value <- as.character(study_value)
if (anyNA(study_value) || any(!nzchar(study_value))) {
stop("`study` cannot contain missing or empty labels.", call. = FALSE)
}
study_value <- make.unique(study_value)
effect_value <- .r4vn_meta_eval(expr$effect, data, env, "effect")
getv <- function(nm) .r4vn_meta_eval(expr[[nm]], data, env, nm)
values <- lapply(
c("se", "lower", "upper", "a", "b", "c", "d",
"event1", "n1", "event0", "n0",
"mean1", "sd1", "mean0", "sd0", "event", "n", "time",
"time1", "time0"),
getv
)
names(values) <- c(
"se", "lower", "upper", "a", "b", "c", "d",
"event1", "n1", "event0", "n0",
"mean1", "sd1", "mean0", "sd0", "event", "n", "time",
"time1", "time0"
)
values$effect <- effect_value
for (nm in names(values)) {
if (!is.null(values[[nm]]) && length(values[[nm]]) != nrow(data)) {
stop("`", nm, "` must have one value per row of `data`.",
call. = FALSE)
}
}
abcd_present <- vapply(values[c("a", "b", "c", "d")],
Negate(is.null), logical(1))
if (any(abcd_present) && !all(abcd_present)) {
stop("Supply all four binary cells: `a`, `b`, `c`, and `d`.",
call. = FALSE)
}
if (all(abcd_present)) {
old_binary <- vapply(values[c("event1", "n1", "event0", "n0")],
Negate(is.null), logical(1))
if (any(old_binary)) {
stop("Use either `a,b,c,d` or `event1,n1,event0,n0`, not both.",
call. = FALSE)
}
values$event1 <- as.numeric(values$a)
values$n1 <- as.numeric(values$a) + as.numeric(values$b)
values$event0 <- as.numeric(values$c)
values$n0 <- as.numeric(values$c) + as.numeric(values$d)
}
inferred <- if (all(abcd_present) ||
all(vapply(values[c("event1", "n1", "event0", "n0")],
Negate(is.null), logical(1)))) {
"OR"
} else if (all(vapply(values[c("event1", "time1", "event0", "time0")],
Negate(is.null), logical(1)))) {
"IRR"
} else if (all(vapply(values[c("mean1", "sd1", "n1", "mean0", "sd0", "n0")],
Negate(is.null), logical(1)))) {
"MD"
} else if (all(vapply(values[c("event", "n")],
Negate(is.null), logical(1)))) {
"PROP"
} else if (all(vapply(values[c("event", "time")],
Negate(is.null), logical(1)))) {
"RATE"
} else NULL
type <- .r4vn_meta_effect_type(
or, rr, rd, hr, irr, md, smd, prop, rate, cor,
!is.null(effect_value), inferred = inferred
)
config <- .r4vn_meta_config(type, transform)
if (identical(config$scale, "ft")) {
warning(
paste(
"Freeman-Tukey transformed proportions are supported, but",
"`transform=\"logit\"` is the R4VN default and is generally preferred."
),
call. = FALSE
)
}
by_names <- unique(c(
.r4vn_meta_varnames(expr$by, data, env, "by", TRUE),
.r4vn_meta_varnames(expr$subgroup, data, env, "subgroup", TRUE)
))
by_name <- if (length(by_names)) by_names[1L] else NULL
cumulative_name <- .r4vn_meta_varname(
expr$cumulative, data, env, "cumulative", TRUE
)
reg_names <- unique(c(.r4vn_meta_reg_names(reg, data),
.r4vn_meta_reg_names(moderator, data)))
raw_required <- switch(
type,
OR = c("event1", "n1", "event0", "n0"),
RR = c("event1", "n1", "event0", "n0"),
RD = c("event1", "n1", "event0", "n0"),
IRR = c("event1", "time1", "event0", "time0"),
MD = c("mean1", "sd1", "n1", "mean0", "sd0", "n0"),
SMD = c("mean1", "sd1", "n1", "mean0", "sd0", "n0"),
PROP = c("event", "n"),
RATE = c("event", "time"),
COR = c("effect", "n"),
HR = character(),
GENERIC = character()
)
use_generic <- !is.null(effect_value) &&
(type %in% c("HR", "OR", "RR", "IRR", "RD", "MD", "SMD", "GENERIC") ||
(type == "COR" && is.null(values$n)) ||
type %in% c("PROP", "RATE"))
if (!use_generic) {
absent <- raw_required[
vapply(raw_required, function(z) is.null(values[[z]]), logical(1))
]
if (length(absent)) {
stop("Missing required input: ", paste(absent, collapse = ", "), ".",
call. = FALSE)
}
}
if (use_generic) {
yi <- .r4vn_meta_transform(as.numeric(effect_value), config)
if (!is.null(values$se)) {
sei <- as.numeric(values$se)
if (any(sei <= 0, na.rm = TRUE)) {
stop("`se` must be greater than 0.", call. = FALSE)
}
} else if (!is.null(values$lower) && !is.null(values$upper)) {
lower_t <- .r4vn_meta_transform(as.numeric(values$lower), config)
upper_t <- .r4vn_meta_transform(as.numeric(values$upper), config)
if (any(lower_t >= upper_t, na.rm = TRUE)) {
stop("Each `lower` value must be below its `upper` value.",
call. = FALSE)
}
zcrit <- stats::qnorm(1 - (1 - ci) / 2)
sei <- (upper_t - lower_t) / (2 * zcrit)
} else {
stop("Generic `effect` input requires `se` or both `lower` and `upper`.",
call. = FALSE)
}
vi <- sei^2
} else {
es <- .r4vn_meta_make_es(type, config, values, cc, zero)
if (is.null(es)) {
stop("This input combination is not supported.", call. = FALSE)
}
yi <- as.numeric(es$yi)
vi <- as.numeric(es$vi)
}
analysis <- data.frame(
.row = seq_len(nrow(data)),
.study = study_value,
yi = yi,
vi = vi,
stringsAsFactors = FALSE
)
if (type %in% c("PROP", "COR") && !is.null(values$n)) {
analysis$.n_for_back <- as.numeric(values$n)
} else {
analysis$.n_for_back <- rep(NA_real_, nrow(data))
}
extra_names <- unique(c(by_names, cumulative_name, reg_names))
extra_names <- extra_names[!is.na(extra_names) & nzchar(extra_names)]
if (length(extra_names)) {
for (nm in extra_names) analysis[[nm]] <- data[[nm]]
}
if (!is.null(values$event1)) analysis$.event1 <- as.numeric(values$event1)
if (!is.null(values$n1)) analysis$.n1 <- as.numeric(values$n1)
if (!is.null(values$event0)) analysis$.event0 <- as.numeric(values$event0)
if (!is.null(values$n0)) analysis$.n0 <- as.numeric(values$n0)
if (!is.null(values$mean1)) analysis$.mean1 <- as.numeric(values$mean1)
if (!is.null(values$sd1)) analysis$.sd1 <- as.numeric(values$sd1)
if (!is.null(values$mean0)) analysis$.mean0 <- as.numeric(values$mean0)
if (!is.null(values$sd0)) analysis$.sd0 <- as.numeric(values$sd0)
if (!is.null(values$event)) analysis$.event <- as.numeric(values$event)
if (!is.null(values$n)) analysis$.n <- as.numeric(values$n)
if (!is.null(values$time)) analysis$.time <- as.numeric(values$time)
if (!is.null(values$time1)) analysis$.time1 <- as.numeric(values$time1)
if (!is.null(values$time0)) analysis$.time0 <- as.numeric(values$time0)
if (all(abcd_present)) {
analysis$.a <- as.numeric(values$a)
analysis$.b <- as.numeric(values$b)
analysis$.c <- as.numeric(values$c)
analysis$.d <- as.numeric(values$d)
}
keep <- is.finite(analysis$yi) & is.finite(analysis$vi) & analysis$vi > 0
if (sum(keep) < 2L) {
stop(
"At least two studies with finite effect estimates and variances are required.",
call. = FALSE
)
}
excluded <- analysis[!keep, c(".study", "yi", "vi"), drop = FALSE]
analysis <- analysis[keep, , drop = FALSE]
rownames(analysis) <- NULL
# Base vector subsetting can drop custom variable/value-label attributes,
# especially when users attach `label` and `labels` without a labelled
# vector class. Rebuild analysis variables from the source data and restore
# these attributes after study exclusion so subgroup output keeps labels.
if (length(extra_names)) {
for (nm in extra_names) {
source_variable <- data[[nm]]
filtered_variable <- source_variable[keep]
for (attribute_name in c("label", "labels", "value.labels")) {
attribute_value <- attr(source_variable, attribute_name, exact = TRUE)
if (!is.null(attribute_value)) {
attr(filtered_variable, attribute_name) <- attribute_value
}
}
analysis[[nm]] <- filtered_variable
}
}
variable_labels <- if (length(extra_names)) {
stats::setNames(
vapply(
extra_names,
function(nm) .r4vn_meta_variable_label(analysis[[nm]], nm),
character(1)
),
extra_names
)
} else character()
test_method <- .r4vn_meta_test_method(small, nrow(analysis), random = random)
random_model <- if (random) {
.r4vn_meta_model(analysis$yi, analysis$vi, method, ci = ci,
fixed = FALSE, test = small)
} else NULL
fixed_model <- if (fixed) {
.r4vn_meta_model(analysis$yi, analysis$vi, method, ci = ci,
fixed = TRUE, test = "z")
} else NULL
primary <- if (!is.null(random_model)) random_model else fixed_model
pred <- .r4vn_meta_prediction(primary, config, ci, analysis$.n_for_back)
overall <- list(
estimate = pred$estimate,
lower = pred$lower,
upper = pred$upper,
pi_lower = pred$pi_lower,
pi_upper = pred$pi_upper,
p = as.numeric(primary$pval)[1L],
statistic = if (!is.null(primary$zval)) {
as.numeric(primary$zval)[1L]
} else NA_real_,
df = if (test_method %in% c("t", "knha", "adhoc")) {
as.numeric((primary$ddf %||% (primary$k - primary$p))[1L])
} else NA_real_,
significant = is.finite(as.numeric(primary$pval)[1L]) &&
as.numeric(primary$pval)[1L] < (1 - ci)
)
subgroup_all <- .r4vn_meta_subgroups_all(
analysis, by_names, method, small, ci, config, digit, p_digit,
fixed = !random, include_interpretation = isTRUE(interpretation)
)
subgroup <- subgroup_all$primary
regression <- .r4vn_meta_regression(
analysis, reg_names, method, small, ci, digit, p_digit, config,
fixed = !random
)
do_prediction <- isTRUE(prediction)
do_bias <- isTRUE(bias)
do_leaveout <- isTRUE(leaveout)
do_influence <- isTRUE(influence)
bias_result <- if (do_bias) {
.r4vn_meta_bias(
primary, primary$k, p_digit, bias_methods, config, ci,
analysis$.n_for_back, strict = strict
)
} else NULL
leave_result <- if (do_leaveout) {
.r4vn_meta_leaveout(primary, config, digit, analysis$.n_for_back)
} else NULL
influence_result <- if (do_influence) {
.r4vn_meta_influence(primary, analysis$.study)
} else NULL
cumulative_result <- .r4vn_meta_cumulative(
analysis, cumulative_name, method, small, ci, config, digit,
fixed = !random
)
small_sample <- if (isTRUE(random)) {
.r4vn_meta_small_sample(
analysis, method, ci, config, analysis$.n_for_back,
requested = small, selected = test_method, digits = digit
)
} else NULL
primary_subgroup <- if (!is.null(by_name)) {
.r4vn_meta_group_values(analysis[[by_name]])
} else NULL
primary_subgroup_title <- if (!is.null(by_name)) {
.r4vn_meta_variable_label(analysis[[by_name]], by_name)
} else "Subgroup"
main_table <- .r4vn_meta_study_table(
analysis, primary, config, digit,
subgroup = primary_subgroup,
subgroup_title = primary_subgroup_title,
ci = ci
)
pooled_rows <- list()
add_pool <- function(label, model) {
pp <- .r4vn_meta_prediction(model, config, ci, analysis$.n_for_back)
row <- as.list(rep("", ncol(main_table)))
names(row) <- names(main_table)
if ("Study" %in% names(row)) row$Study <- label
row[["Effect (95% CI)"]] <- .r4vn_meta_ci_text(
pp$estimate, pp$lower, pp$upper, digit
)
row$Weight <- "100.0%"
as.data.frame(row, stringsAsFactors = FALSE, check.names = FALSE)
}
if (fixed) {
pooled_rows[[length(pooled_rows) + 1L]] <-
add_pool("Overall (fixed)", fixed_model)
}
if (random) {
pooled_rows[[length(pooled_rows) + 1L]] <-
add_pool("Overall (random)", random_model)
}
pool <- do.call(rbind, pooled_rows)
main_table <- rbind(main_table, pool)
rownames(main_table) <- NULL
heterogeneity <- .r4vn_meta_heterogeneity(primary, p_digit)
significance_tests <- .r4vn_meta_significance_tests(
primary, overall, test_method, ci, p_digit,
include_interpretation = isTRUE(interpretation)
)
if (do_prediction &&
is.finite(overall$pi_lower) && is.finite(overall$pi_upper)) {
heterogeneity <- rbind(
heterogeneity,
data.frame(
Statistic = paste0(round(ci * 100), "% prediction interval"),
Value = paste0(
.r4vn_meta_fmt(overall$pi_lower, digit), "-",
.r4vn_meta_fmt(overall$pi_upper, digit)
),
stringsAsFactors = FALSE
)
)
}
if (is.null(title) || !length(title) || is.na(title[1L]) ||
!nzchar(as.character(title)[1L])) {
title <- paste0("Meta-analysis of ", config$label)
} else {
title <- as.character(title)[1L]
}
report_dir <- tempfile("r4vn-meta-")
dir.create(report_dir, recursive = TRUE, showWarnings = FALSE)
binary_table <- .r4vn_meta_binary_table(analysis)
overview <- data.frame(
Studies = primary$k,
Effect = config$label,
Model = if (random) "Random effects" else "Fixed effect",
`Tau-squared method` = if (random) method else "Not applicable",
Inference = test_method,
stringsAsFactors = FALSE, check.names = FALSE
)
out <- list(
data = main_table,
table = main_table,
tables = list(
Overview = overview, Main = main_table,
Statistical_tests = significance_tests,
Heterogeneity = heterogeneity
),
title = title,
config = config,
effect_type = type,
analysis_data = analysis,
excluded = excluded,
model = primary,
primary_model = primary,
random_model = random_model,
fixed_model = fixed_model,
random = random,
fixed = fixed,
method = method,
hk = test_method %in% c("knha", "adhoc"),
small = small,
test_method = test_method,
small_sample = small_sample,
prediction = do_prediction,
overall = overall,
significance_tests = significance_tests,
heterogeneity = heterogeneity,
subgroup = subgroup,
subgroups = subgroup_all$analyses,
subgroup_tables = subgroup_all$wide,
subgroup_long = subgroup_all$table,
regression = regression,
bias = bias_result,
leaveout = leave_result,
influence = influence_result,
cumulative = cumulative_result,
by = by_name,
by_vars = by_names,
reg = reg_names,
cumulative_var = cumulative_name,
full = full,
profile = profile,
report = report_text,
interpretation = NULL,
overview = overview,
ci = ci,
digits = digit,
p_digits = p_digit,
plots = list(),
plot_titles = list(),
plot_display = plot_display,
plot_args = plot_args,
report_dir = report_dir,
call = match.call()
)
if (!is.null(binary_table)) out$tables$Binary_2x2 <- binary_table
if (length(subgroup_all$wide)) {
out$tables$Subgroup <- subgroup_all$wide[[1L]]
if (length(subgroup_all$wide) > 1L) {
additional_subgroups <- subgroup_all$wide[-1L]
for (nm in names(additional_subgroups)) {
out$tables[[paste0("Subgroup_", nm)]] <- additional_subgroups[[nm]]
}
}
}
if (!is.null(subgroup_all$test)) out$tables$Subgroup_test <- subgroup_all$test
if (!is.null(regression$table)) {
out$tables$Meta_regression <- regression$table
}
if (!is.null(regression$univariable)) {
out$tables$Moderator_univariable <- regression$univariable
}
if (!is.null(regression$statistics)) {
out$tables$Meta_regression_statistics <- regression$statistics
}
if (!is.null(bias_result$table)) {
out$tables$Publication_bias <- bias_result$table
}
if (!is.null(bias_result$adjusted)) {
out$tables$Publication_bias_adjusted <- bias_result$adjusted
}
if (!is.null(small_sample$comparison)) {
out$tables$Small_sample_inference <- small_sample$comparison
}
if (!is.null(leave_result$table)) {
out$tables$Leave_one_out <- leave_result$table
}
if (!is.null(influence_result$table)) {
out$tables$Influence <- influence_result$table
}
if (!is.null(cumulative_result$table)) {
out$tables$Cumulative <- cumulative_result$table
}
out$estimates <- list(
overall = overall,
subgroup = subgroup_all$table,
meta_regression = if (is.null(regression)) NULL else regression$table,
publication_bias_adjusted = if (is.null(bias_result)) NULL else bias_result$adjusted,
cumulative = if (is.null(cumulative_result)) NULL else cumulative_result$table
)
out$tests <- list(
overall = significance_tests[1L, , drop = FALSE],
heterogeneity = significance_tests[2L, , drop = FALSE],
significance = significance_tests,
subgroup = subgroup_all$test,
moderators = if (is.null(regression)) NULL else regression$statistics,
small_study_effects = if (is.null(bias_result)) NULL else bias_result$table
)
out$diagnostics <- list(
small_sample = small_sample,
publication_bias = bias_result,
leave_one_out = leave_result,
influence = influence_result
)
out$models <- list(
primary = primary,
random = random_model,
fixed = fixed_model,
subgroups = subgroup_all$analyses,
meta_regression = if (is.null(regression)) NULL else regression$model,
trimfill = if (is.null(bias_result)) NULL else bias_result$trimfill,
selection = if (is.null(bias_result)) NULL else bias_result$selection
)
out$metadata <- list(
profile = profile, effect_type = type, effect_label = config$label,
studies = primary$k, excluded = nrow(excluded), by = by_names,
moderators = reg_names, method = method, inference = test_method,
ci = ci, continuity_correction = cc, zero = zero,
variable_labels = variable_labels,
plot_display = plot_display,
input = if (all(abcd_present)) "a,b,c,d" else if (use_generic) "generic" else "raw"
)
class(out) <- c("r4vn_meta", "r4vn_tab")
out$results_text <- .r4vn_meta_results_text(out)
if (isTRUE(report_text)) {
out$interpretation <- .r4vn_meta_interpretation(out)
out$tables$Interpretation <- out$interpretation
}
if (isTRUE(plot)) {
plot_opts <- function(type) {
nested <- any(names(plot_args) %in% c(
"all", "forest", "funnel", "trimfill", "leaveout", "influence",
"baujat", "radial", "labbe", "cumulative", "bubble"
))
opts <- if (nested) {
utils::modifyList(plot_args$all %||% list(), plot_args[[type]] %||% list())
} else plot_args
opts[c("x", "type", "file")] <- NULL
opts
}
save_plot <- function(type, path, extra = list(), force = list(),
object = out) {
opts <- utils::modifyList(extra, plot_opts(type))
opts <- utils::modifyList(opts, force)
args <- c(list(x = object, type = type, file = path), opts)
ok <- tryCatch(
{
do.call(.r4vn_meta_plot_file, args)
TRUE
},
error = function(e) {
if (isTRUE(strict)) stop(e)
warning("Meta-analysis plot `", type, "` was unavailable: ",
conditionMessage(e), call. = FALSE)
FALSE
}
)
if (!isTRUE(ok) && file.exists(path)) unlink(path)
invisible(ok)
}
subgroup_plot_object <- function(subgroup_analysis, group_name) {
idx <- subgroup_analysis$indices[[group_name]]
if (is.null(idx) || sum(idx) < 2L) return(NULL)
subgroup_object <- out
subgroup_object$analysis_data <- analysis[idx, , drop = FALSE]
subgroup_object$primary_model <- subgroup_analysis$fits[[group_name]]
subgroup_object$model <- subgroup_object$primary_model
subgroup_object$random_model <- if (random) {
subgroup_object$primary_model
} else NULL
subgroup_object$fixed_model <- if (fixed && !random) {
subgroup_object$primary_model
} else NULL
subgroup_object$by <- NULL
subgroup_object$title <- paste0(
"Subgroup forest plot: ",
subgroup_analysis$label %||% subgroup_analysis$variable,
" = ", group_name
)
subgroup_object
}
out$plots$forest <- file.path(report_dir, "forest.png")
out$plot_titles$forest <- "Forest plot"
save_plot("forest", out$plots$forest)
if (length(out$subgroups)) {
subgroup_number <- 0L
for (moderator_name in names(out$subgroups)) {
subgroup_analysis <- out$subgroups[[moderator_name]]
if (!length(subgroup_analysis$fits)) next
for (group_name in names(subgroup_analysis$fits)) {
subgroup_number <- subgroup_number + 1L
subgroup_object <- subgroup_plot_object(
subgroup_analysis, group_name
)
if (is.null(subgroup_object)) next
subgroup_title <- subgroup_object$title
plot_name <- paste0("forest_subgroup_", subgroup_number)
out$plots[[plot_name]] <- file.path(
report_dir, paste0(plot_name, ".png")
)
out$plot_titles[[plot_name]] <- subgroup_title
save_plot(
"forest", out$plots[[plot_name]],
force = list(title = subgroup_title), object = subgroup_object
)
}
}
}
if (do_bias && primary$k >= 3L) {
out$plots$funnel <- file.path(report_dir, "funnel.png")
out$plot_titles$funnel <- "Contour-enhanced funnel plot"
save_plot("funnel", out$plots$funnel, list(contour = TRUE))
if (!is.null(bias_result$trimfill)) {
out$plots$trimfill <- file.path(report_dir, "trimfill.png")
out$plot_titles$trimfill <- "Trim-and-fill funnel plot"
save_plot("trimfill", out$plots$trimfill, list(contour = TRUE))
}
}
if (do_leaveout && !is.null(leave_result$table)) {
out$plots$leaveout <- file.path(report_dir, "leaveout.png")
out$plot_titles$leaveout <- "Leave-one-out analysis"
save_plot("leaveout", out$plots$leaveout)
}
if (do_influence && !is.null(influence_result$table)) {
out$plots$influence <- file.path(report_dir, "influence.png")
out$plot_titles$influence <- "Influence diagnostics"
save_plot("influence", out$plots$influence)
out$plots$baujat <- file.path(report_dir, "baujat.png")
out$plot_titles$baujat <- "Baujat plot"
save_plot("baujat", out$plots$baujat)
}
if (!is.null(cumulative_result$table)) {
out$plots$cumulative <- file.path(report_dir, "cumulative.png")
out$plot_titles$cumulative <- paste0(
"Cumulative meta-analysis: ",
cumulative_result$label %||% cumulative_name
)
save_plot(
"cumulative", out$plots$cumulative,
list(title = out$plot_titles$cumulative)
)
}
if (!is.null(regression$model) && length(reg_names) == 1L) {
moderator_label <- .r4vn_meta_variable_label(
analysis[[reg_names[1L]]], reg_names[1L]
)
out$plots$bubble <- file.path(report_dir, "bubble.png")
out$plot_titles$bubble <- paste0("Meta-regression: ", moderator_label)
save_plot("bubble", out$plots$bubble,
list(moderator = reg_names[1L],
title = out$plot_titles$bubble))
}
if (identical(profile, "full")) {
out$plots$radial <- file.path(report_dir, "radial.png")
out$plot_titles$radial <- "Radial plot"
save_plot("radial", out$plots$radial)
if (!is.null(binary_table)) {
out$plots$labbe <- file.path(report_dir, "labbe.png")
out$plot_titles$labbe <- "L'Abbe plot"
save_plot("labbe", out$plots$labbe)
}
}
available <- vapply(out$plots, file.exists, logical(1))
out$plots <- out$plots[available]
out$plot_titles <- out$plot_titles[names(out$plots)]
}
out$html <- .r4vn_meta_report_html(out)
out$table_html <- .r4vn_meta_html_table(out$table, out$title)
out$file <- file.path(report_dir, "index.html")
writeLines(enc2utf8(out$html), out$file, useBytes = TRUE)
if (isTRUE(plot) && interactive() && length(plot_display)) {
display_types <- if ("all" %in% plot_display) {
intersect(names(out$plots), plot_types)
} else plot_display
missing_display <- setdiff(display_types, names(out$plots))
if (length(missing_display)) {
warning(
"Plot-pane figure(s) unavailable for this analysis: ",
paste(missing_display, collapse = ", "), ".",
call. = FALSE
)
}
display_types <- intersect(display_types, names(out$plots))
for (display_type in display_types) {
display_args <- c(
list(x = out, type = display_type), plot_opts(display_type)
)
display_args[c("file", "width", "height", "res")] <- NULL
tryCatch(
do.call(.r4vn_meta_draw, display_args),
error = function(e) {
if (isTRUE(strict)) stop(e)
warning(
"Meta-analysis Plot-pane figure `", display_type,
"` was unavailable: ", conditionMessage(e), call. = FALSE
)
}
)
if (identical(display_type, "forest") && length(out$subgroups)) {
subgroup_drawn <- 0L
for (moderator_name in names(out$subgroups)) {
subgroup_analysis <- out$subgroups[[moderator_name]]
if (!length(subgroup_analysis$fits)) next
for (group_name in names(subgroup_analysis$fits)) {
subgroup_object <- subgroup_plot_object(
subgroup_analysis, group_name
)
if (is.null(subgroup_object)) next
subgroup_args <- c(
list(x = subgroup_object, type = "forest"),
plot_opts("forest")
)
subgroup_args[c("file", "width", "height", "res")] <- NULL
# A subgroup title must remain explicit even if a common forest
# title was supplied through plot_args.
subgroup_args$title <- subgroup_object$title
drawn <- tryCatch(
{
do.call(.r4vn_meta_draw, subgroup_args)
TRUE
},
error = function(e) {
if (isTRUE(strict)) stop(e)
warning(
"Subgroup forest plot `", subgroup_object$title,
"` was unavailable in the Plot pane: ",
conditionMessage(e), call. = FALSE
)
FALSE
}
)
if (isTRUE(drawn)) subgroup_drawn <- subgroup_drawn + 1L
}
}
if (subgroup_drawn > 0L) {
message(
"R4VN added the overall forest plot and ", subgroup_drawn,
" subgroup forest plot(s) to Plot history. Use Previous/Next ",
"to review them."
)
}
}
}
}
if (!is.null(export) || !is.null(file)) {
export_format <- export
if (is.null(export_format)) {
export_format <- tools::file_ext(as.character(file)[1L])
if (!nzchar(export_format)) export_format <- "html"
}
out$export <- .r4vn_export_meta(
out, export = export_format, file = file,
open = open, title = title
)
}
if (isTRUE(show)) .r4vn_meta_open(out$file)
invisible(out)
}
#' Print or reopen an R4VN meta-analysis
#'
#' @param x An object created by `tabmeta()`.
#' @param ... Additional arguments ignored.
#' @return The input object invisibly.
#' @export
print.r4vn_meta <- function(x, ...) {
if (!inherits(x, "r4vn_meta")) {
stop("`x` must be created by `tabmeta()`.", call. = FALSE)
}
if (!is.null(x$file) && file.exists(x$file)) {
.r4vn_meta_open(x$file)
} else {
cat(x$title, "\n")
print(x$table, row.names = FALSE)
cat("\n")
print(x$significance_tests, row.names = FALSE)
cat("\n")
print(x$heterogeneity, row.names = FALSE)
}
invisible(x)
}
#' Convert an R4VN meta-analysis to a data frame
#'
#' @param x An object created by `tabmeta()`.
#' @param row.names Ignored.
#' @param optional Ignored.
#' @param ... Additional arguments ignored.
#' @return The main publication table.
#' @export
as.data.frame.r4vn_meta <- function(x, row.names = NULL,
optional = FALSE, ...) {
x$table
}
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 .r4vn_meta_subgroups .r4vn_meta_metric_ci_text .r4vn_meta_heterogeneity_ci .r4vn_meta_heterogeneity .r4vn_meta_binary_table .r4vn_meta_study_table .r4vn_meta_term_labels .r4vn_meta_group_values .r4vn_meta_group_spec .r4vn_meta_variable_label .r4vn_meta_weights .r4vn_meta_prediction .r4vn_meta_model .r4vn_meta_make_es .r4vn_meta_ci_text_compact .r4vn_meta_fmt_compact .r4vn_meta_ci_text .r4vn_meta_p .r4vn_meta_fmt .r4vn_meta_back_fun .r4vn_meta_back .r4vn_meta_transform .r4vn_meta_config .r4vn_meta_effect_type .r4vn_meta_test_method .r4vn_meta_profile .r4vn_meta_varnames .r4vn_meta_reg_names .r4vn_meta_varname .r4vn_meta_eval .r4vn_meta_escape .r4vn_meta_num1 .r4vn_meta_flag .r4vn_meta_need
# ============================================================================
# R4VN tabmeta: publication-ready meta-analysis
# Backend: metafor (optional Suggests dependency)
# ============================================================================
.r4vn_meta_need <- function(package = "metafor") {
if (!requireNamespace(package, quietly = TRUE)) {
stop(
"Package `", package, "` is required for this analysis. Install it with ",
"install.packages(\"", package, "\").",
call. = FALSE
)
}
invisible(TRUE)
}
.r4vn_meta_flag <- function(x, name) {
if (!is.logical(x) || length(x) != 1L || is.na(x)) {
stop("`", name, "` must be TRUE or FALSE.", call. = FALSE)
}
x
}
.r4vn_meta_num1 <- function(x, name, lower = -Inf, upper = Inf,
inclusive_lower = TRUE, inclusive_upper = TRUE) {
if (!is.numeric(x) || length(x) != 1L || is.na(x) || !is.finite(x)) {
stop("`", name, "` must be one finite numeric value.", call. = FALSE)
}
ok_lower <- if (inclusive_lower) x >= lower else x > lower
ok_upper <- if (inclusive_upper) x <= upper else x < upper
if (!ok_lower || !ok_upper) {
stop("`", name, "` is outside the allowed range.", call. = FALSE)
}
x
}
.r4vn_meta_escape <- function(x) {
x <- as.character(x)
x <- gsub("&", "&", x, fixed = TRUE)
x <- gsub("<", "<", x, fixed = TRUE)
x <- gsub(">", ">", x, fixed = TRUE)
x <- gsub('"', """, x, fixed = TRUE)
gsub("'", "'", x, fixed = TRUE)
}
.r4vn_meta_eval <- function(expr, data, env, arg, required = FALSE) {
if (identical(expr, quote(NULL)) || is.null(expr)) {
if (required) stop("`", arg, "` is required.", call. = FALSE)
return(NULL)
}
value <- tryCatch(
eval(expr, envir = data, enclos = env),
error = function(e) tryCatch(eval(expr, envir = env),
error = function(e2) NULL)
)
if (is.null(value) && required) {
stop("Could not evaluate `", arg, "`.", call. = FALSE)
}
if (!is.null(value) && length(value) == 1L && nrow(data) > 1L &&
is.character(value) && value %in% names(data)) {
value <- data[[value]]
}
value
}
.r4vn_meta_varname <- function(expr, data, env, arg, allow_null = TRUE) {
if (identical(expr, quote(NULL)) || is.null(expr)) {
if (allow_null) return(NULL)
stop("`", arg, "` is required.", call. = FALSE)
}
if (is.symbol(expr)) {
nm <- as.character(expr)
if (nm %in% names(data)) return(nm)
}
value <- tryCatch(eval(expr, envir = data, enclos = env),
error = function(e) NULL)
if (is.character(value) && length(value) == 1L && value %in% names(data)) {
return(value)
}
stop("`", arg, "` must identify one variable in `data`.", call. = FALSE)
}
.r4vn_meta_reg_names <- function(reg, data) {
if (is.null(reg)) return(character())
if (inherits(reg, "r4vn_vars")) {
out <- as.character(reg$variable)
} else if (is.character(reg)) {
out <- as.character(reg)
} else {
stop("`reg` must be created using `vars()` or be a character vector.",
call. = FALSE)
}
out <- unique(out[nzchar(out)])
absent <- setdiff(out, names(data))
if (length(absent)) {
stop("Variables in `reg` not found in `data`: ",
paste(absent, collapse = ", "), ".", call. = FALSE)
}
out
}
.r4vn_meta_varnames <- function(expr, data, env, arg, allow_null = TRUE) {
if (identical(expr, quote(NULL)) || is.null(expr)) {
if (allow_null) return(character())
stop("`", arg, "` is required.", call. = FALSE)
}
if (is.symbol(expr)) {
nm <- as.character(expr)
if (nm %in% names(data)) return(nm)
}
value <- tryCatch(eval(expr, envir = data, enclos = env),
error = function(e) NULL)
if (inherits(value, "r4vn_vars")) {
out <- as.character(value$variable)
} else if (is.character(value)) {
out <- as.character(value)
} else {
stop("`", arg, "` must identify variables in `data` or use `vars()`.",
call. = FALSE)
}
out <- unique(out[!is.na(out) & nzchar(out)])
absent <- setdiff(out, names(data))
if (length(absent)) {
stop("Variables in `", arg, "` not found in `data`: ",
paste(absent, collapse = ", "), ".", call. = FALSE)
}
out
}
.r4vn_meta_profile <- function(profile, full, prediction, bias, leaveout,
influence, plot) {
profile <- match.arg(profile, c("auto", "brief", "full", "custom"))
if (!is.null(full)) {
.r4vn_meta_flag(full, "full")
profile <- if (isTRUE(full)) "full" else "custom"
}
defaults <- switch(
profile,
auto = list(prediction = TRUE, bias = TRUE, leaveout = TRUE,
influence = TRUE, plot = TRUE),
brief = list(prediction = FALSE, bias = FALSE, leaveout = FALSE,
influence = FALSE, plot = TRUE),
full = list(prediction = TRUE, bias = TRUE, leaveout = TRUE,
influence = TRUE, plot = TRUE),
custom = list(prediction = FALSE, bias = FALSE, leaveout = FALSE,
influence = FALSE, plot = FALSE)
)
supplied <- list(prediction = prediction, bias = bias, leaveout = leaveout,
influence = influence, plot = plot)
for (nm in names(defaults)) {
if (!is.null(supplied[[nm]])) {
.r4vn_meta_flag(supplied[[nm]], nm)
defaults[[nm]] <- supplied[[nm]]
}
}
c(defaults, list(profile = profile, full = identical(profile, "full")))
}
.r4vn_meta_test_method <- function(small, k, random = TRUE) {
small <- match.arg(tolower(as.character(small)[1L]),
c("auto", "adhoc", "knha", "t", "z"))
if (!isTRUE(random)) return("z")
if (identical(small, "auto")) {
if (k <= 10L) "adhoc" else "z"
} else small
}
.r4vn_meta_effect_type <- function(or, rr, rd, hr, irr, md, smd, prop,
rate, cor, effect_present,
inferred = NULL) {
flags <- c(OR = or, RR = rr, RD = rd, HR = hr, IRR = irr, MD = md,
SMD = smd, PROP = prop, RATE = rate, COR = cor)
selected <- names(flags)[flags]
if (length(selected) > 1L) {
stop("Choose only one effect option: `or`, `rr`, `rd`, `hr`, `irr`, ",
"`md`, `smd`, `prop`, `rate`, or `cor`.", call. = FALSE)
}
if (!length(selected)) {
if (effect_present) return("GENERIC")
if (!is.null(inferred) && length(inferred)) return(inferred[1L])
stop("Choose an effect type such as `or=TRUE`, `rr=TRUE`, `hr=TRUE`, ",
"`md=TRUE`, `smd=TRUE`, `prop=TRUE`, `rate=TRUE`, or `cor=TRUE`.",
call. = FALSE)
}
selected
}
.r4vn_meta_config <- function(type, transform = NULL) {
if (type %in% c("OR", "RR", "HR", "IRR")) {
return(list(measure = type, ref = 1, label = type, scale = "log"))
}
if (type == "RD") {
return(list(measure = "RD", ref = 0, label = "Risk difference",
scale = "identity"))
}
if (type == "MD") {
return(list(measure = "MD", ref = 0, label = "Mean difference",
scale = "identity"))
}
if (type == "SMD") {
return(list(measure = "SMD", ref = 0,
label = "Standardized mean difference", scale = "identity"))
}
if (type == "COR") {
return(list(measure = "ZCOR", ref = 0, label = "Correlation",
scale = "zcor"))
}
if (type == "RATE") {
return(list(measure = "IRLN", ref = 0, label = "Incidence rate",
scale = "log"))
}
if (type == "PROP") {
if (is.null(transform)) transform <- "logit"
transform <- tolower(as.character(transform)[1L])
transform <- switch(
transform,
logit = "logit",
arcsine = "arcsine",
asin = "arcsine",
ft = "ft",
`freeman-tukey` = "ft",
none = "none",
stop("`transform` for proportions must be `logit`, `arcsine`, `ft`, or `none`.",
call. = FALSE)
)
measure <- switch(transform, logit = "PLO", arcsine = "PAS",
ft = "PFT", none = "PR")
return(list(measure = measure, ref = 0, label = "Proportion",
scale = transform))
}
list(measure = "GEN", ref = 0, label = "Effect", scale = "identity")
}
.r4vn_meta_transform <- function(x, config) {
if (config$scale == "log") {
if (any(x <= 0, na.rm = TRUE)) {
stop("Ratio/rate estimates must be greater than 0.", call. = FALSE)
}
return(log(x))
}
if (config$scale == "zcor") {
if (any(abs(x) >= 1, na.rm = TRUE)) {
stop("Correlation estimates must be strictly between -1 and 1.",
call. = FALSE)
}
return(atanh(x))
}
x
}
.r4vn_meta_back <- function(x, config, ni = NULL) {
if (config$scale == "log") return(exp(x))
if (config$scale == "zcor") return(tanh(x))
if (config$scale == "logit") return(stats::plogis(x))
if (config$scale == "arcsine") return(sin(x)^2)
if (config$scale == "ft") {
.r4vn_meta_need()
if (is.null(ni)) return(rep(NA_real_, length(x)))
if (length(x) == length(ni) && length(x) > 1L) {
return(mapply(
function(z, nn) metafor::transf.ipft(z, targs = list(ni = nn)),
x, ni, USE.NAMES = FALSE
))
}
return(metafor::transf.ipft.hm(x, targs = list(ni = ni)))
}
x
}
.r4vn_meta_back_fun <- function(x) {
config <- x$config
ni <- x$analysis_data$.n_for_back
if (config$scale == "log") return(exp)
if (config$scale == "zcor") return(tanh)
if (config$scale == "logit") return(stats::plogis)
if (config$scale == "arcsine") return(function(z) sin(z)^2)
if (config$scale == "ft") {
return(function(z) metafor::transf.ipft.hm(z, targs = list(ni = ni)))
}
identity
}
.r4vn_meta_fmt <- function(x, digits = 2) {
if (!length(x)) return(character(0))
out <- rep("", length(x))
ok <- !is.na(x) & is.finite(x)
if (any(ok)) {
out[ok] <- formatC(
x[ok], format = "f", digits = digits, big.mark = ","
)
}
out
}
.r4vn_meta_p <- function(x, digits = 3) {
if (!length(x) || is.na(x) || !is.finite(x)) return("")
lim <- 10^(-digits)
if (x < lim) paste0("<", formatC(lim, format = "f", digits = digits)) else
formatC(x, format = "f", digits = digits)
}
.r4vn_meta_ci_text <- function(est, lo, hi, digits = 2) {
if (any(!is.finite(c(est, lo, hi)))) return("")
paste0(.r4vn_meta_fmt(est, digits), " (",
.r4vn_meta_fmt(lo, digits), "-", .r4vn_meta_fmt(hi, digits), ")")
}
.r4vn_meta_fmt_compact <- function(x, digits = 2) {
if (!length(x)) return(character())
out <- .r4vn_meta_fmt(x, digits)
finite <- is.finite(x)
scientific <- finite & x != 0 &
(abs(x) >= 1e6 | abs(x) < 10^(-digits))
if (any(scientific)) {
out[scientific] <- formatC(
x[scientific], format = "e", digits = digits
)
}
out
}
.r4vn_meta_ci_text_compact <- function(est, lo, hi, digits = 2) {
if (any(!is.finite(c(est, lo, hi)))) return("")
paste0(
.r4vn_meta_fmt_compact(est, digits), " (",
.r4vn_meta_fmt_compact(lo, digits), "-",
.r4vn_meta_fmt_compact(hi, digits), ")"
)
}
.r4vn_meta_make_es <- function(type, config, values, cc, zero) {
.r4vn_meta_need()
drop00 <- identical(zero, "exclude")
if (type %in% c("OR", "RR", "RD")) {
ai <- values$event1
ci0 <- values$event0
n1 <- values$n1
n0 <- values$n0
if (any(ai < 0 | ci0 < 0 | n1 <= 0 | n0 <= 0, na.rm = TRUE) ||
any(ai > n1 | ci0 > n0, na.rm = TRUE)) {
stop("Binary counts must satisfy 0 <= event <= n and n > 0.",
call. = FALSE)
}
return(metafor::escalc(
measure = config$measure,
ai = ai, bi = n1 - ai,
ci = ci0, di = n0 - ci0,
add = cc, to = if (cc == 0) "none" else "only0",
drop00 = drop00
))
}
if (type %in% c("MD", "SMD")) {
if (any(values$n1 <= 1 | values$n0 <= 1, na.rm = TRUE) ||
any(values$sd1 < 0 | values$sd0 < 0, na.rm = TRUE)) {
stop("Continuous input requires n1,n0 > 1 and non-negative SDs.",
call. = FALSE)
}
return(metafor::escalc(
measure = config$measure,
m1i = values$mean1, sd1i = values$sd1, n1i = values$n1,
m2i = values$mean0, sd2i = values$sd0, n2i = values$n0
))
}
if (type == "PROP") {
if (any(values$event < 0 | values$n <= 0 | values$event > values$n,
na.rm = TRUE)) {
stop("Proportion input requires 0 <= event <= n and n > 0.",
call. = FALSE)
}
return(metafor::escalc(
measure = config$measure, xi = values$event, ni = values$n,
add = cc, to = if (cc == 0) "none" else "only0"
))
}
if (type == "RATE") {
if (any(values$event < 0 | values$time <= 0, na.rm = TRUE)) {
stop("Rate input requires event >= 0 and person-time > 0.",
call. = FALSE)
}
return(metafor::escalc(
measure = "IRLN", xi = values$event, ti = values$time,
add = cc, to = if (cc == 0) "none" else "only0"
))
}
if (type == "IRR") {
if (any(values$event1 < 0 | values$event0 < 0 |
values$time1 <= 0 | values$time0 <= 0, na.rm = TRUE)) {
stop("IRR input requires non-negative events and positive person-time.",
call. = FALSE)
}
return(metafor::escalc(
measure = "IRR",
x1i = values$event1, x2i = values$event0,
t1i = values$time1, t2i = values$time0,
add = cc, to = if (cc == 0) "none" else "only0"
))
}
if (type == "COR") {
if (any(abs(values$effect) >= 1, na.rm = TRUE) ||
any(values$n <= 3, na.rm = TRUE)) {
stop("Correlation input requires -1 < effect < 1 and n > 3.",
call. = FALSE)
}
return(metafor::escalc(
measure = "ZCOR", ri = values$effect, ni = values$n
))
}
NULL
}
.r4vn_meta_model <- function(yi, vi, method, hk = NULL, ci, fixed = FALSE,
mods = NULL, test = NULL) {
.r4vn_meta_need()
if (is.null(test)) {
test <- if (!fixed && isTRUE(hk)) "knha" else "z"
}
test <- .r4vn_meta_test_method(test, length(yi), random = !fixed)
args <- list(
yi = yi,
vi = vi,
method = if (fixed) "FE" else method,
test = test,
level = ci * 100
)
if (!is.null(mods)) args$mods <- mods
do.call(metafor::rma.uni, args)
}
.r4vn_meta_prediction <- function(model, config, ci, ni = NULL) {
pred <- stats::predict(model, level = ci * 100)
list(
estimate = .r4vn_meta_back(as.numeric(pred$pred), config, ni),
lower = .r4vn_meta_back(as.numeric(pred$ci.lb), config, ni),
upper = .r4vn_meta_back(as.numeric(pred$ci.ub), config, ni),
pi_lower = if (!is.null(pred$pi.lb))
.r4vn_meta_back(as.numeric(pred$pi.lb), config, ni) else NA_real_,
pi_upper = if (!is.null(pred$pi.ub))
.r4vn_meta_back(as.numeric(pred$pi.ub), config, ni) else NA_real_
)
}
.r4vn_meta_weights <- function(model) {
w <- tryCatch(stats::weights(model), error = function(e) NULL)
if (is.null(w)) return(rep(NA_real_, model$k))
as.numeric(w)
}
.r4vn_meta_variable_label <- function(x, fallback) {
label <- attr(x, "label", exact = TRUE)
if (is.null(label) || !length(label) || is.na(label[1L]) ||
!nzchar(trimws(as.character(label)[1L]))) {
return(as.character(fallback)[1L])
}
as.character(label)[1L]
}
.r4vn_meta_group_spec <- function(x) {
observed <- x[!is.na(x)]
raw <- if (is.factor(x)) {
levels(droplevels(x))
} else {
unique(as.character(observed))
}
raw <- raw[raw %in% as.character(observed)]
display <- raw
value_labels <- attr(x, "labels", exact = TRUE)
if (is.null(value_labels)) {
value_labels <- attr(x, "value.labels", exact = TRUE)
}
if (!is.null(value_labels) && length(value_labels) &&
!is.null(names(value_labels))) {
labelled_values <- as.character(unname(value_labels))
matched <- match(raw, labelled_values)
use <- !is.na(matched) & nzchar(names(value_labels)[matched])
display[use] <- names(value_labels)[matched[use]]
}
duplicate <- duplicated(display) | duplicated(display, fromLast = TRUE)
display[duplicate] <- paste0(display[duplicate], " (", raw[duplicate], ")")
list(raw = raw, display = display)
}
.r4vn_meta_group_values <- function(x) {
spec <- .r4vn_meta_group_spec(x)
out <- as.character(x)
matched <- match(out, spec$raw)
use <- !is.na(matched)
out[use] <- spec$display[matched[use]]
out[is.na(x)] <- NA_character_
out
}
.r4vn_meta_term_labels <- function(terms, data, variables) {
out <- as.character(terms)
normalized <- tolower(gsub("[()[:space:]]", "", out))
out[normalized %in% c("intrcpt", "intercept")] <- "Intercept"
variables <- unique(as.character(variables))
variables <- variables[variables %in% names(data)]
if (!length(variables)) return(out)
variables <- variables[order(nchar(variables), decreasing = TRUE)]
for (nm in variables) {
variable_label <- .r4vn_meta_variable_label(data[[nm]], nm)
exact <- terms == nm
out[exact] <- variable_label
prefixed <- !exact & startsWith(as.character(terms), nm)
if (any(prefixed)) {
suffix <- substring(as.character(terms)[prefixed], nchar(nm) + 1L)
suffix <- sub("^[.:_]+", "", suffix)
out[prefixed] <- ifelse(
nzchar(suffix), paste0(variable_label, ": ", suffix), variable_label
)
}
}
out
}
.r4vn_meta_study_table <- function(d, model, config, digits,
subgroup = NULL, subgroup_title = "Subgroup",
ci = 0.95) {
back <- .r4vn_meta_back_fun(list(config = config, analysis_data = d))
zcrit <- stats::qnorm(1 - (1 - ci) / 2)
est <- back(d$yi)
lo <- back(d$yi - zcrit * sqrt(d$vi))
hi <- back(d$yi + zcrit * sqrt(d$vi))
wt <- .r4vn_meta_weights(model)
if (length(wt) != nrow(d)) wt <- rep(NA_real_, nrow(d))
out <- data.frame(
Study = d$.study,
stringsAsFactors = FALSE, check.names = FALSE
)
if (all(c(".event1", ".n1", ".event0", ".n0") %in% names(d))) {
out[["Group 1"]] <- paste0(
.r4vn_meta_fmt(d$.event1, 0), "/", .r4vn_meta_fmt(d$.n1, 0)
)
out[["Group 0"]] <- paste0(
.r4vn_meta_fmt(d$.event0, 0), "/", .r4vn_meta_fmt(d$.n0, 0)
)
} else if (all(c(".event1", ".time1", ".event0", ".time0") %in% names(d))) {
out[["Group 1"]] <- paste0(
.r4vn_meta_fmt(d$.event1, 0), "/",
.r4vn_meta_fmt(d$.time1, digits)
)
out[["Group 0"]] <- paste0(
.r4vn_meta_fmt(d$.event0, 0), "/",
.r4vn_meta_fmt(d$.time0, digits)
)
} else if (all(c(".mean1", ".sd1", ".n1", ".mean0", ".sd0", ".n0") %in%
names(d))) {
out[["Group 1"]] <- paste0(
.r4vn_meta_fmt(d$.mean1, digits), " (",
.r4vn_meta_fmt(d$.sd1, digits), "); n=",
.r4vn_meta_fmt(d$.n1, 0)
)
out[["Group 0"]] <- paste0(
.r4vn_meta_fmt(d$.mean0, digits), " (",
.r4vn_meta_fmt(d$.sd0, digits), "); n=",
.r4vn_meta_fmt(d$.n0, 0)
)
} else if (all(c(".event", ".n") %in% names(d))) {
out[["Events/Total"]] <- paste0(
.r4vn_meta_fmt(d$.event, 0), "/", .r4vn_meta_fmt(d$.n, 0)
)
} else if (all(c(".event", ".time") %in% names(d))) {
out[["Events/Person-time"]] <- paste0(
.r4vn_meta_fmt(d$.event, 0), "/",
.r4vn_meta_fmt(d$.time, digits)
)
}
out[["Effect (95% CI)"]] <- mapply(
.r4vn_meta_ci_text, est, lo, hi,
MoreArgs = list(digits = digits), USE.NAMES = FALSE
)
out[["Weight"]] <- ifelse(
is.finite(wt), paste0(.r4vn_meta_fmt(wt, 1), "%"), ""
)
if (!is.null(subgroup)) {
subgroup_data <- data.frame(subgroup, stringsAsFactors = FALSE,
check.names = FALSE)
names(subgroup_data) <- as.character(subgroup_title)[1L]
out <- cbind(subgroup_data, out)
}
out
}
.r4vn_meta_binary_table <- function(d) {
required <- c(".event1", ".n1", ".event0", ".n0")
if (!all(required %in% names(d))) return(NULL)
a <- d$.event1
b <- d$.n1 - d$.event1
c0 <- d$.event0
d0 <- d$.n0 - d$.event0
data.frame(
Study = d$.study,
`a: event, group 1` = a,
`b: non-event, group 1` = b,
`c: event, group 0` = c0,
`d: non-event, group 0` = d0,
`Group 1 total` = d$.n1,
`Group 0 total` = d$.n0,
stringsAsFactors = FALSE, check.names = FALSE
)
}
.r4vn_meta_heterogeneity <- function(model, p_digits = 3) {
data.frame(
Statistic = c("Studies", "Q", "df", "p heterogeneity", "I-squared",
"H-squared", "tau-squared", "tau"),
Value = c(
as.character(model$k),
.r4vn_meta_fmt(model$QE, 2),
as.character(model$k - model$p),
.r4vn_meta_p(model$QEp, p_digits),
paste0(.r4vn_meta_fmt(model$I2, 1), "%"),
.r4vn_meta_fmt(model$H2, 2),
.r4vn_meta_fmt(model$tau2, 4),
.r4vn_meta_fmt(sqrt(model$tau2), 4)
),
stringsAsFactors = FALSE, check.names = FALSE
)
}
.r4vn_meta_heterogeneity_ci <- function(model, ci = 0.95) {
empty <- list(
i2 = c(estimate = NA_real_, lower = NA_real_, upper = NA_real_),
tau2 = c(estimate = NA_real_, lower = NA_real_, upper = NA_real_)
)
result <- tryCatch(
suppressWarnings(stats::confint(model, level = ci * 100)),
error = function(e) NULL
)
if (is.null(result)) return(empty)
random <- if (is.list(result) && !is.null(result$random)) {
result$random
} else if (is.matrix(result) || is.data.frame(result)) {
result
} else NULL
if (is.null(random)) return(empty)
random <- as.data.frame(random, check.names = FALSE)
if (!nrow(random) || ncol(random) < 3L) return(empty)
row_key <- tolower(gsub("[^a-zA-Z0-9]", "", rownames(random)))
column_key <- tolower(gsub("[^a-zA-Z0-9]", "", names(random)))
estimate_column <- match("estimate", column_key)
lower_column <- match("cilb", column_key)
upper_column <- match("ciub", column_key)
if (anyNA(c(estimate_column, lower_column, upper_column))) {
estimate_column <- 1L
lower_column <- 2L
upper_column <- 3L
}
extract <- function(key) {
row <- match(key, row_key)
if (is.na(row)) return(c(estimate = NA_real_, lower = NA_real_,
upper = NA_real_))
values <- suppressWarnings(as.numeric(random[
row, c(estimate_column, lower_column, upper_column), drop = TRUE
]))
stats::setNames(values, c("estimate", "lower", "upper"))
}
list(i2 = extract("i2"), tau2 = extract("tau2"))
}
.r4vn_meta_metric_ci_text <- function(estimate, interval, digits = 2,
suffix = "") {
if (!is.finite(estimate)) return("")
if (length(interval) >= 3L &&
all(is.finite(interval[c("lower", "upper")]))) {
return(paste0(
.r4vn_meta_fmt(estimate, digits), suffix, " (",
.r4vn_meta_fmt(interval[["lower"]], digits), suffix, "-",
.r4vn_meta_fmt(interval[["upper"]], digits), suffix, ")"
))
}
paste0(.r4vn_meta_fmt(estimate, digits), suffix)
}
.r4vn_meta_subgroups <- function(d, by_name, method, test, ci, config,
digits, p_digits, fixed = FALSE,
include_interpretation = FALSE) {
if (is.null(by_name)) return(NULL)
g <- d[[by_name]]
group_spec <- .r4vn_meta_group_spec(g)
lev <- group_spec$raw
display_lev <- group_spec$display
variable_label <- .r4vn_meta_variable_label(g, by_name)
if (length(lev) < 2L) return(NULL)
fits <- list()
indices <- list()
rows <- list()
retained_raw <- retained_display <- character()
for (i in seq_along(lev)) {
z <- lev[i]
display_z <- display_lev[i]
idx <- !is.na(g) & as.character(g) == z
if (sum(idx) < 2L) next
fit <- .r4vn_meta_model(d$yi[idx], d$vi[idx], method, ci = ci,
test = test, fixed = fixed)
fits[[display_z]] <- fit
indices[[display_z]] <- idx
retained_raw <- c(retained_raw, z)
retained_display <- c(retained_display, display_z)
pr <- .r4vn_meta_prediction(fit, config, ci, d$.n_for_back[idx])
heterogeneity_ci <- .r4vn_meta_heterogeneity_ci(fit, ci)
fit_test <- tolower(as.character(fit$test %||% test)[1L])
test_symbol <- if (fit_test %in% c("knha", "adhoc", "t")) "t" else "z"
effect_statistic <- as.numeric(fit$zval)[1L]
effect_df <- if (test_symbol == "t") {
as.numeric((fit$ddf %||% (fit$k - fit$p))[1L])
} else NA_real_
effect_p <- as.numeric(fit$pval)[1L]
prediction_interval <- if (!fixed && is.finite(pr$pi_lower) &&
is.finite(pr$pi_upper)) {
paste0(
.r4vn_meta_fmt(pr$pi_lower, digits), "-",
.r4vn_meta_fmt(pr$pi_upper, digits)
)
} else ""
subgroup_row <- data.frame(
Subgroup = display_z,
Studies = fit$k,
`Effect (95% CI)` = .r4vn_meta_ci_text(
pr$estimate, pr$lower, pr$upper, digits
),
`Test statistic` = if (is.finite(effect_statistic)) {
paste0(test_symbol, " = ", .r4vn_meta_fmt(effect_statistic, 2))
} else "",
`df effect` = if (is.finite(effect_df)) {
.r4vn_meta_fmt(effect_df, 0)
} else "",
`p-value effect` = .r4vn_meta_p(effect_p, p_digits),
`Prediction interval` = prediction_interval,
`Q heterogeneity` = .r4vn_meta_fmt(fit$QE, 2),
`df heterogeneity` = as.character(max(0, fit$k - fit$p)),
`p-value heterogeneity` = .r4vn_meta_p(fit$QEp, p_digits),
`I-squared (95% CI)` = .r4vn_meta_metric_ci_text(
fit$I2, heterogeneity_ci$i2, digits = 1, suffix = "%"
),
`Tau-squared (95% CI)` = .r4vn_meta_metric_ci_text(
fit$tau2, heterogeneity_ci$tau2, digits = 4
),
stringsAsFactors = FALSE, check.names = FALSE
)
confidence_percent <- round(ci * 100)
names(subgroup_row)[names(subgroup_row) == "Effect (95% CI)"] <-
paste0("Effect (", confidence_percent, "% CI)")
names(subgroup_row)[names(subgroup_row) == "Prediction interval"] <-
paste0(confidence_percent, "% prediction interval")
names(subgroup_row)[names(subgroup_row) == "I-squared (95% CI)"] <-
paste0("I-squared (", confidence_percent, "% CI)")
names(subgroup_row)[names(subgroup_row) == "Tau-squared (95% CI)"] <-
paste0("Tau-squared (", confidence_percent, "% CI)")
if (isTRUE(include_interpretation)) {
effect_significant <- is.finite(effect_p) && effect_p < (1 - ci)
heterogeneity_significant <- is.finite(fit$QEp) &&
fit$QEp < (1 - ci)
effect_conclusion <- if (!is.finite(effect_p)) {
"the pooled-effect significance test was unavailable"
} else if (effect_significant) {
"the pooled effect was statistically significant"
} else {
"the pooled effect was not statistically significant"
}
heterogeneity_conclusion <- if (!is.finite(fit$QEp)) {
"the heterogeneity significance test was unavailable."
} else if (heterogeneity_significant) {
"heterogeneity was statistically significant."
} else {
"heterogeneity was not statistically significant."
}
subgroup_row[["Statistically significant"]] <- if (
is.finite(effect_p)
) if (effect_significant) "Yes" else "No" else ""
subgroup_row[["Conclusion"]] <- paste0(
toupper(substr(effect_conclusion, 1L, 1L)),
substring(effect_conclusion, 2L), "; ",
heterogeneity_conclusion
)
}
rows[[display_z]] <- subgroup_row
}
table <- if (length(rows)) do.call(rbind, rows) else NULL
if (!is.null(table)) rownames(table) <- NULL
complete_group <- !is.na(g)
moddata <- data.frame(
.g = factor(as.character(g[complete_group]))
)
mm <- stats::model.matrix(~ .g, data = moddata)
testfit <- if (ncol(mm) > 1L) tryCatch(
.r4vn_meta_model(
d$yi[complete_group], d$vi[complete_group], method, ci = ci,
mods = mm[, -1L, drop = FALSE], test = test, fixed = fixed
),
error = function(e) NULL
) else NULL
test <- if (is.null(testfit)) NULL else data.frame(
Statistic = c("Q between subgroups", "df", "p"),
Value = c(
.r4vn_meta_fmt(testfit$QM, 2),
as.character(testfit$m),
.r4vn_meta_p(testfit$QMp, p_digits)
),
stringsAsFactors = FALSE, check.names = FALSE
)
list(
fits = fits, indices = indices, table = table, test = test,
test_model = testfit, levels = retained_display,
raw_levels = retained_raw, variable = by_name, label = variable_label
)
}
.r4vn_meta_subgroup_transpose <- function(table) {
if (is.null(table) || !is.data.frame(table) || !nrow(table) ||
!"Subgroup" %in% names(table)) return(NULL)
metrics <- setdiff(names(table), "Subgroup")
out <- data.frame(
Statistic = metrics, stringsAsFactors = FALSE, check.names = FALSE
)
for (i in seq_len(nrow(table))) {
group_name <- as.character(table$Subgroup[i])
values <- vapply(
table[i, metrics, drop = FALSE],
function(value) {
if (!length(value) || is.na(value[1L])) "" else as.character(value[1L])
},
character(1)
)
out[[group_name]] <- unname(values)
}
out
}
.r4vn_meta_significance_tests <- function(model, overall, test_method,
ci = 0.95, p_digits = 3,
include_interpretation = FALSE) {
alpha <- 1 - ci
pooled_stat <- if (!is.null(model$zval) && length(model$zval)) {
as.numeric(model$zval[1L])
} else {
NA_real_
}
pooled_df <- if (test_method %in% c("t", "knha", "adhoc")) {
candidate <- model$ddf %||% (model$k - model$p)
as.numeric(candidate[1L])
} else {
NA_real_
}
pooled_p <- as.numeric(overall$p)[1L]
heterogeneity_p <- as.numeric(model$QEp)[1L]
pooled_sig <- is.finite(pooled_p) && pooled_p < alpha
heterogeneity_sig <- is.finite(heterogeneity_p) &&
heterogeneity_p < alpha
pooled_symbol <- if (test_method %in% c("t", "knha", "adhoc")) {
"t"
} else {
"z"
}
out <- data.frame(
Test = c("Pooled effect", "Heterogeneity (Cochran's Q)"),
`Test statistic` = c(
if (is.finite(pooled_stat)) {
paste0(pooled_symbol, " = ", .r4vn_meta_fmt(pooled_stat, 2))
} else "",
if (is.finite(model$QE)) {
paste0("Q = ", .r4vn_meta_fmt(model$QE, 2))
} else ""
),
df = c(
if (is.finite(pooled_df)) .r4vn_meta_fmt(pooled_df, 0) else "",
.r4vn_meta_fmt(max(0, model$k - model$p), 0)
),
`p-value` = c(
.r4vn_meta_p(pooled_p, p_digits),
.r4vn_meta_p(heterogeneity_p, p_digits)
),
stringsAsFactors = FALSE, check.names = FALSE
)
if (!isTRUE(include_interpretation)) return(out)
out[["Statistically significant"]] <- c(
if (is.finite(pooled_p)) if (pooled_sig) "Yes" else "No" else "",
if (is.finite(heterogeneity_p)) {
if (heterogeneity_sig) "Yes" else "No"
} else ""
)
out[["Conclusion"]] <- c(
if (!is.finite(pooled_p)) {
"Statistical significance of the pooled estimate was not available."
} else if (pooled_sig) {
"The pooled estimate is statistically significant."
} else {
"The pooled estimate is not statistically significant."
},
if (!is.finite(heterogeneity_p)) {
"Statistical significance of heterogeneity was not available."
} else if (heterogeneity_sig) {
"There is statistically significant between-study heterogeneity."
} else {
"There is no statistically significant evidence of between-study heterogeneity."
}
)
out
}
.r4vn_meta_subgroups_all <- function(d, by_names, method, test, ci, config,
digits, p_digits, fixed = FALSE,
include_interpretation = FALSE) {
if (!length(by_names)) return(list(primary = NULL, analyses = list(),
wide = list(), table = NULL,
test = NULL))
analyses <- lapply(by_names, function(nm) {
.r4vn_meta_subgroups(
d, nm, method, test, ci, config, digits, p_digits, fixed = fixed,
include_interpretation = include_interpretation
)
})
names(analyses) <- by_names
analyses <- Filter(Negate(is.null), analyses)
wide <- lapply(analyses, function(result) {
.r4vn_meta_subgroup_transpose(result$table)
})
wide <- Filter(Negate(is.null), wide)
bind_component <- function(component) {
rows <- lapply(names(analyses), function(nm) {
z <- analyses[[nm]][[component]]
if (is.null(z) || !is.data.frame(z) || !nrow(z)) return(NULL)
data.frame(Moderator = analyses[[nm]]$label %||% nm, z,
stringsAsFactors = FALSE,
check.names = FALSE)
})
rows <- Filter(Negate(is.null), rows)
if (!length(rows)) NULL else {
out <- do.call(rbind, rows)
rownames(out) <- NULL
out
}
}
list(
primary = if (length(analyses)) analyses[[1L]] else NULL,
analyses = analyses,
wide = wide,
table = bind_component("table"),
test = bind_component("test")
)
}
.r4vn_meta_regression <- function(d, reg_names, method, test, ci, digits,
p_digits, config = NULL, fixed = FALSE) {
if (!length(reg_names)) return(NULL)
md <- d[, reg_names, drop = FALSE]
keep <- stats::complete.cases(md)
if (sum(keep) < max(3L, length(reg_names) + 2L)) {
return(list(status = "Not enough complete studies for meta-regression."))
}
mm <- stats::model.matrix(~ ., data = md[keep, , drop = FALSE])
if (ncol(mm) < 2L) {
return(list(status = "No estimable moderator terms."))
}
fit <- .r4vn_meta_model(
d$yi[keep], d$vi[keep], method, ci = ci,
mods = mm[, -1L, drop = FALSE], test = test, fixed = fixed
)
sm <- summary(fit)
beta <- as.numeric(sm$beta)
se <- as.numeric(sm$se)
p <- as.numeric(sm$pval)
rn <- rownames(sm$beta)
if (is.null(rn)) rn <- c("Intercept", colnames(mm)[-1L])
display_rn <- .r4vn_meta_term_labels(rn, d, reg_names)
lower_beta <- if (!is.null(sm$ci.lb)) as.numeric(sm$ci.lb) else {
beta - stats::qnorm(1 - (1 - ci) / 2) * se
}
upper_beta <- if (!is.null(sm$ci.ub)) as.numeric(sm$ci.ub) else {
beta + stats::qnorm(1 - (1 - ci) / 2) * se
}
tab <- data.frame(
Moderator = display_rn,
`Beta (95% CI)` = mapply(
.r4vn_meta_ci_text_compact,
beta, lower_beta, upper_beta,
MoreArgs = list(digits = digits), USE.NAMES = FALSE
),
`p-value` = vapply(p, .r4vn_meta_p, character(1), digits = p_digits),
stringsAsFactors = FALSE, check.names = FALSE
)
confidence_percent <- round(ci * 100)
names(tab)[names(tab) == "Beta (95% CI)"] <- paste0(
"Beta (", confidence_percent, "% CI)"
)
if (!is.null(config) && identical(config$scale, "log")) {
tab[["Ratio of effects (95% CI)"]] <- mapply(
.r4vn_meta_ci_text_compact,
exp(beta), exp(lower_beta), exp(upper_beta),
MoreArgs = list(digits = digits), USE.NAMES = FALSE
)
names(tab)[names(tab) == "Ratio of effects (95% CI)"] <- paste0(
"Ratio of effects (", confidence_percent, "% CI)"
)
}
stats_tab <- data.frame(
Statistic = c("QM", "df", "p moderators", "Residual I-squared",
"Residual tau-squared"),
Value = c(
.r4vn_meta_fmt(fit$QM, 2),
as.character(fit$m),
.r4vn_meta_p(fit$QMp, p_digits),
paste0(.r4vn_meta_fmt(fit$I2, 1), "%"),
.r4vn_meta_fmt(fit$tau2, 4)
),
stringsAsFactors = FALSE, check.names = FALSE
)
univariable <- lapply(reg_names, function(nm) {
mdi <- d[, nm, drop = FALSE]
keep_i <- stats::complete.cases(mdi)
if (sum(keep_i) < 3L) return(NULL)
mm_i <- stats::model.matrix(~ ., data = mdi[keep_i, , drop = FALSE])
if (ncol(mm_i) < 2L) return(NULL)
fit_i <- tryCatch(
.r4vn_meta_model(
d$yi[keep_i], d$vi[keep_i], method, ci = ci,
mods = mm_i[, -1L, drop = FALSE], test = test, fixed = fixed
),
error = function(e) NULL
)
if (is.null(fit_i)) return(NULL)
data.frame(
Moderator = .r4vn_meta_variable_label(d[[nm]], nm),
Studies = sum(keep_i),
QM = .r4vn_meta_fmt(fit_i$QM, 2),
df = fit_i$m,
`p-value` = .r4vn_meta_p(fit_i$QMp, p_digits),
`Residual I-squared` = paste0(.r4vn_meta_fmt(fit_i$I2, 1), "%"),
stringsAsFactors = FALSE, check.names = FALSE
)
})
univariable <- Filter(Negate(is.null), univariable)
list(model = fit, table = tab, statistics = stats_tab,
univariable = if (length(univariable)) do.call(rbind, univariable) else NULL,
keep = keep, model_matrix = mm)
}
.r4vn_meta_bias_methods <- function(methods, profile = "auto") {
if (is.null(methods) || !length(methods) || identical(methods, "auto")) {
methods <- if (identical(profile, "full")) {
c("egger", "begg", "trimfill", "failsafe", "selection")
} else c("egger", "begg", "trimfill")
}
methods <- unique(tolower(as.character(methods)))
aliases <- c(`rank` = "begg", `trim-and-fill` = "trimfill",
`fail-safe` = "failsafe", `selmodel` = "selection")
hit <- methods %in% names(aliases)
methods[hit] <- unname(aliases[methods[hit]])
allowed <- c("egger", "begg", "trimfill", "failsafe", "selection")
bad <- setdiff(methods, allowed)
if (length(bad)) {
stop("Unsupported publication-bias method: ", paste(bad, collapse = ", "),
".", call. = FALSE)
}
methods
}
.r4vn_meta_bias <- function(model, k, p_digits, methods, config, ci,
ni = NULL, strict = FALSE) {
rows <- list()
adjusted <- list()
notes <- character()
add_test <- function(test, statistic = "", p = "", result = "Available") {
rows[[length(rows) + 1L]] <<- data.frame(
Test = test, Statistic = statistic, `p-value` = p, Result = result,
stringsAsFactors = FALSE, check.names = FALSE
)
}
safe <- function(expr, label) {
warning_messages <- character()
value <- tryCatch(
withCallingHandlers(
expr,
warning = function(w) {
warning_messages <<- c(warning_messages, conditionMessage(w))
invokeRestart("muffleWarning")
}
),
error = function(e) {
if (isTRUE(strict)) stop(e)
notes <<- c(
notes,
paste0(label, " was unavailable: ", conditionMessage(e))
)
NULL
}
)
if (length(warning_messages)) {
notes <<- c(
notes,
paste0(
label, " warning: ",
paste(unique(warning_messages), collapse = " ")
)
)
}
value
}
egger <- begg <- trimfill <- failsafe <- selection <- NULL
if ("egger" %in% methods) {
if (k < 10L) {
add_test("Egger regression test", result = "Not performed (<10 studies)")
} else {
egger <- safe(metafor::regtest(model, model = "lm", predictor = "sei"),
"Egger test")
if (!is.null(egger)) add_test(
"Egger regression test", .r4vn_meta_fmt(egger$zval, 3),
.r4vn_meta_p(egger$pval, p_digits)
)
}
}
if ("begg" %in% methods) {
if (k < 10L) {
add_test("Begg-Mazumdar rank correlation",
result = "Not performed (<10 studies)")
} else {
begg <- safe(metafor::ranktest(model), "Rank-correlation test")
if (!is.null(begg)) {
stat <- if (!is.null(begg$tau)) begg$tau else begg$zval
add_test("Begg-Mazumdar rank correlation",
.r4vn_meta_fmt(stat, 3),
.r4vn_meta_p(begg$pval, p_digits))
}
}
}
if ("trimfill" %in% methods && k >= 3L) {
trimfill <- safe(metafor::trimfill(model), "Trim-and-fill analysis")
if (!is.null(trimfill)) {
pr <- .r4vn_meta_prediction(trimfill, config, ci, ni)
adjusted[[length(adjusted) + 1L]] <- data.frame(
Method = "Trim-and-fill",
`Adjusted effect (95% CI)` = .r4vn_meta_ci_text(
pr$estimate, pr$lower, pr$upper, 3
),
`Imputed studies` = as.integer(trimfill$k0 %||% 0L),
stringsAsFactors = FALSE, check.names = FALSE
)
}
}
if ("failsafe" %in% methods && k >= 3L) {
failsafe <- safe(
metafor::fsn(x = model$yi, vi = model$vi, type = "Rosenberg"),
"Fail-safe N"
)
if (!is.null(failsafe)) add_test(
"Rosenberg fail-safe N",
as.character(failsafe$fsnum %||% ""),
.r4vn_meta_p(failsafe$pval %||% NA_real_, p_digits),
"Exploratory file-drawer analysis"
)
}
if ("selection" %in% methods && k >= 10L) {
selection <- safe(
metafor::selmodel(model, type = "stepfun", steps = c(.025, .05, .10, .50)),
"Selection model"
)
if (!is.null(selection)) {
b <- as.numeric(selection$beta[1L])
se <- as.numeric(selection$se[1L])
crit <- stats::qnorm(1 - (1 - ci) / 2)
adjusted[[length(adjusted) + 1L]] <- data.frame(
Method = "Step-function selection model",
`Adjusted effect (95% CI)` = .r4vn_meta_ci_text(
.r4vn_meta_back(b, config, ni),
.r4vn_meta_back(b - crit * se, config, ni),
.r4vn_meta_back(b + crit * se, config, ni), 3
),
`Imputed studies` = NA_integer_,
stringsAsFactors = FALSE, check.names = FALSE
)
}
}
if ("selection" %in% methods && k < 10L) {
notes <- c(notes, "Selection models were not fitted with fewer than 10 studies.")
}
table <- if (length(rows)) do.call(rbind, rows) else NULL
adjusted_table <- if (length(adjusted)) do.call(rbind, adjusted) else NULL
list(
status = if (length(notes)) paste(unique(notes), collapse = " ") else NULL,
table = table, adjusted = adjusted_table, methods = methods,
egger = egger, begg = begg, trimfill = trimfill,
failsafe = failsafe, selection = selection
)
}
.r4vn_meta_small_sample <- function(d, method, ci, config, ni, requested,
selected, digits = 2) {
methods <- unique(c("z", "knha", "adhoc", selected))
rows <- lapply(methods, function(test) {
fit <- tryCatch(
.r4vn_meta_model(d$yi, d$vi, method, ci = ci, test = test),
error = function(e) NULL
)
if (is.null(fit)) return(NULL)
pr <- .r4vn_meta_prediction(fit, config, ci, ni)
data.frame(
Inference = switch(test, z = "Normal approximation",
t = "t distribution",
knha = "Hartung-Knapp",
adhoc = "Modified Hartung-Knapp (adhoc)"),
Selected = ifelse(identical(test, selected), "Yes", "No"),
df = if (test == "z") NA_real_ else max(1, fit$k - fit$p),
`Effect (95% CI)` = .r4vn_meta_ci_text(
pr$estimate, pr$lower, pr$upper, digits
),
stringsAsFactors = FALSE, check.names = FALSE
)
})
rows <- Filter(Negate(is.null), rows)
comparison <- if (length(rows)) do.call(rbind, rows) else NULL
status <- if (nrow(d) <= 10L && identical(selected, "adhoc")) {
paste(
"Few-study inference is active. The modified Hartung-Knapp method",
"prevents the adjustment from producing smaller standard errors than",
"the conventional random-effects analysis."
)
} else if (nrow(d) <= 10L) {
paste0(
"Only ", nrow(d), " studies are available; inference method `",
selected, "` was selected. The sensitivity comparison is shown."
)
} else if (identical(requested, "auto")) {
"The normal approximation is used because more than 10 studies are available."
} else {
paste0("Inference method `", selected,
"` was explicitly selected; the sensitivity comparison is shown.")
}
list(
studies = nrow(d), requested = requested, selected = selected,
threshold = 10L,
status = status,
comparison = comparison
)
}
.r4vn_meta_leaveout <- function(model, config, digits, ni = NULL) {
if (model$k < 3L) return(NULL)
z <- tryCatch(metafor::leave1out(model), error = function(e) NULL)
if (is.null(z)) return(NULL)
zz <- as.data.frame(z)
if (!all(c("estimate", "ci.lb", "ci.ub") %in% names(zz))) return(NULL)
est <- .r4vn_meta_back(zz$estimate, config, ni)
lo <- .r4vn_meta_back(zz$ci.lb, config, ni)
hi <- .r4vn_meta_back(zz$ci.ub, config, ni)
slab <- rownames(zz)
if (is.null(slab) || !length(slab)) {
slab <- paste0("Study ", seq_len(nrow(zz)))
}
i2 <- if ("I2" %in% names(zz)) zz$I2 else rep(NA_real_, nrow(zz))
out <- data.frame(
`Study omitted` = slab,
`Effect (95% CI)` = mapply(
.r4vn_meta_ci_text, est, lo, hi,
MoreArgs = list(digits = digits), USE.NAMES = FALSE
),
`I-squared` = ifelse(
is.finite(i2), paste0(.r4vn_meta_fmt(i2, 1), "%"), ""
),
stringsAsFactors = FALSE, check.names = FALSE
)
list(raw = z, table = out, estimate = est, lower = lo, upper = hi)
}
.r4vn_meta_influence <- function(model, study) {
if (model$k < 3L) return(NULL)
inf <- tryCatch(stats::influence(model), error = function(e) NULL)
if (is.null(inf)) return(NULL)
d <- inf$inf
if (is.null(d)) return(list(raw = inf, table = NULL))
d <- as.data.frame(d)
pick <- function(candidates) {
hit <- intersect(candidates, names(d))
if (!length(hit)) rep(NA_real_, nrow(d)) else as.numeric(d[[hit[1L]]])
}
cook <- pick(c("cook.d", "cook.d."))
rstudent <- pick(c("rstudent", "rstudent."))
hat <- pick(c("hat", "hat."))
is_infl <- inf$is.infl
if (is.null(is_infl) || length(is_infl) != nrow(d)) {
is_infl <- rep(FALSE, nrow(d))
}
if (length(study) != nrow(d)) study <- rownames(d)
tab <- data.frame(
Study = study,
`Studentized residual` = vapply(
rstudent, .r4vn_meta_fmt, character(1), digits = 3
),
`Cook's distance` = vapply(
cook, .r4vn_meta_fmt, character(1), digits = 3
),
`Hat value` = vapply(
hat, .r4vn_meta_fmt, character(1), digits = 3
),
Influential = ifelse(is_infl, "Yes", "No"),
stringsAsFactors = FALSE, check.names = FALSE
)
list(raw = inf, table = tab, cook = cook, rstudent = rstudent,
hat = hat, influential = is_infl)
}
.r4vn_meta_cumulative <- function(d, order_name, method, test, ci, config,
digits, fixed = FALSE) {
if (is.null(order_name)) return(NULL)
order_label <- .r4vn_meta_variable_label(d[[order_name]], order_name)
ord <- order(d[[order_name]], na.last = NA)
if (length(ord) < 2L) return(NULL)
dd <- d[ord, , drop = FALSE]
rows <- vector("list", nrow(dd) - 1L)
for (i in 2:nrow(dd)) {
fit <- .r4vn_meta_model(
dd$yi[seq_len(i)], dd$vi[seq_len(i)], method, ci = ci, test = test,
fixed = fixed
)
pr <- .r4vn_meta_prediction(
fit, config, ci, dd$.n_for_back[seq_len(i)]
)
rows[[i - 1L]] <- data.frame(
`Up to` = as.character(dd[[order_name]][i]),
Studies = i,
`Effect (95% CI)` = .r4vn_meta_ci_text(
pr$estimate, pr$lower, pr$upper, digits
),
`I-squared` = paste0(.r4vn_meta_fmt(fit$I2, 1), "%"),
Estimate = pr$estimate,
Lower = pr$lower,
Upper = pr$upper,
stringsAsFactors = FALSE, check.names = FALSE
)
}
plot_data <- do.call(rbind, rows)
names(plot_data)[1L] <- paste0("Up to: ", order_label)
tab <- plot_data
tab$Estimate <- tab$Lower <- tab$Upper <- NULL
list(table = tab, plot_data = plot_data, order = order_name,
label = order_label)
}
.r4vn_meta_interpretation <- function(x) {
rows <- list()
add <- function(section, text) {
if (is.null(text) || !length(text) || is.na(text[1L]) ||
!nzchar(trimws(as.character(text)[1L]))) return(invisible(NULL))
rows[[length(rows) + 1L]] <<- data.frame(
Section = section,
Interpretation = as.character(text)[1L],
stringsAsFactors = FALSE, check.names = FALSE
)
invisible(NULL)
}
ptxt <- function(p) {
if (is.null(p) || !length(p) || !is.finite(p[1L])) return("not available")
paste0("p ", if (p[1L] < 10^(-x$p_digits)) {
paste0("< ", formatC(10^(-x$p_digits), format = "f",
digits = x$p_digits))
} else {
paste0("= ", formatC(p[1L], format = "f", digits = x$p_digits))
})
}
ci <- function(est, lo, hi, digits = x$digits) {
paste0(.r4vn_meta_fmt(est, digits), " (", round(x$ci * 100), "% CI ",
.r4vn_meta_fmt(lo, digits), " to ",
.r4vn_meta_fmt(hi, digits), ")")
}
null <- x$config$ref
crosses_null <- function(lo, hi) {
is.finite(lo) && is.finite(hi) && lo <= null && hi >= null
}
pr <- x$overall
h <- x$primary_model
model_name <- if (x$random) "random-effects" else "fixed-effect"
has_comparative_null <- !x$effect_type %in% c("PROP", "RATE")
direction <- if (!is.finite(pr$estimate) || !has_comparative_null) {
""
} else if (pr$estimate < null) {
paste0(" The pooled estimate was below the null value of ", null, ".")
} else if (pr$estimate > null) {
paste0(" The pooled estimate was above the null value of ", null, ".")
} else {
paste0(" The pooled estimate equaled the null value of ", null, ".")
}
null_text <- if (!has_comparative_null) {
""
} else if (crosses_null(pr$lower, pr$upper)) {
" The confidence interval included the null value."
} else {
" The confidence interval excluded the null value."
}
inference_name <- switch(
x$test_method, z = "normal-approximation", t = "t-distribution",
knha = "Hartung-Knapp", adhoc = "modified Hartung-Knapp",
x$test_method
)
add(
"Overall effect",
paste0(
h$k, " studies were included. The ", model_name,
" model estimated a pooled ", x$config$label, " of ",
ci(pr$estimate, pr$lower, pr$upper), "; ", ptxt(pr$p), ". ",
if (is.finite(pr$p) && pr$p < (1 - x$ci)) {
"The pooled estimate was statistically significant."
} else if (is.finite(pr$p)) {
"The pooled estimate was not statistically significant."
} else {
"Statistical significance of the pooled estimate was not available."
},
direction, null_text, " Inference used the ", inference_name,
" method."
)
)
i2 <- as.numeric(h$I2)
i2_description <- if (!is.finite(i2)) {
"not estimable"
} else if (i2 < 25) {
"low"
} else if (i2 < 50) {
"moderate"
} else if (i2 < 75) {
"substantial"
} else {
"considerable"
}
add(
"Heterogeneity",
paste0(
"Between-study heterogeneity was ", i2_description,
" on the I-squared scale (I-squared = ",
.r4vn_meta_fmt(i2, 1), "%; tau-squared = ",
.r4vn_meta_fmt(h$tau2, 4), "). Cochran's Q was ",
.r4vn_meta_fmt(h$QE, 2), " with ", max(0, h$k - h$p),
" degrees of freedom (", ptxt(h$QEp), "). ",
if (is.finite(h$QEp) && h$QEp < (1 - x$ci)) {
"Heterogeneity was statistically significant."
} else if (is.finite(h$QEp)) {
"There was no statistically significant evidence of heterogeneity."
} else {
"Statistical significance of heterogeneity was not available."
}
)
)
if (isTRUE(x$prediction) && isTRUE(x$random) &&
is.finite(pr$pi_lower) && is.finite(pr$pi_upper)) {
pi_text <- if (!has_comparative_null) {
NULL
} else if (crosses_null(pr$pi_lower, pr$pi_upper)) {
"included"
} else {
"excluded"
}
add(
"Prediction interval",
paste0(
"The ", round(x$ci * 100), "% prediction interval was ",
.r4vn_meta_fmt(pr$pi_lower, x$digits), " to ",
.r4vn_meta_fmt(pr$pi_upper, x$digits),
if (is.null(pi_text)) "." else paste0(" and ", pi_text, " the null value."),
" This interval describes the range expected for the ",
"underlying effect in a comparable new study, conditional on the fitted model."
)
)
}
if (length(x$subgroups)) {
for (nm in names(x$subgroups)) {
z <- x$subgroups[[nm]]
fit <- z$test_model
if (is.null(fit)) next
evidence <- if (is.finite(fit$QMp) && fit$QMp < 0.05) {
"There was statistical evidence that pooled effects differed between subgroups."
} else {
"There was no statistical evidence that pooled effects differed between subgroups."
}
add(
paste0("Subgroup analysis: ", z$label %||% nm),
paste0(
evidence, " The test for subgroup differences gave Q = ",
.r4vn_meta_fmt(fit$QM, 2), " with ", fit$m,
" degrees of freedom (", ptxt(fit$QMp), "). This test should be ",
"interpreted cautiously when subgroups contain few studies."
)
)
}
}
if (!is.null(x$regression$model)) {
fit <- x$regression$model
sm <- summary(fit)
beta_names <- rownames(sm$beta)
if (is.null(beta_names)) beta_names <- paste0("coefficient ", seq_along(sm$pval))
display_beta_names <- .r4vn_meta_term_labels(
beta_names, x$analysis_data, x$reg
)
moderator_idx <- !grepl("^(intrcpt|intercept|\\(intercept\\))$",
beta_names, ignore.case = TRUE)
sig <- display_beta_names[
moderator_idx & is.finite(sm$pval) & sm$pval < 0.05
]
sig_text <- if (length(sig)) {
paste0(" Statistically associated coefficient(s): ",
paste(sig, collapse = ", "), ".")
} else {
" No individual moderator coefficient met the 0.05 significance threshold."
}
add(
"Meta-regression",
paste0(
"The joint moderator test gave QM = ", .r4vn_meta_fmt(fit$QM, 2),
" with ", fit$m, " degrees of freedom (", ptxt(fit$QMp), ").",
sig_text, " Residual heterogeneity was I-squared = ",
.r4vn_meta_fmt(fit$I2, 1), "% and tau-squared = ",
.r4vn_meta_fmt(fit$tau2, 4),
". Meta-regression findings are observational across studies and should not be interpreted causally."
)
)
}
if (!is.null(x$small_sample$status)) {
add("Few-study inference", x$small_sample$status)
}
if (!is.null(x$bias)) {
bias_bits <- character()
if (!is.null(x$bias$egger)) {
bias_bits <- c(
bias_bits,
paste0(
"Egger's regression test ",
if (x$bias$egger$pval < 0.05) "detected" else "did not detect",
" funnel asymmetry (", ptxt(x$bias$egger$pval), ")"
)
)
}
if (!is.null(x$bias$begg)) {
bias_bits <- c(
bias_bits,
paste0(
"the rank-correlation test ",
if (x$bias$begg$pval < 0.05) "detected" else "did not detect",
" asymmetry (", ptxt(x$bias$begg$pval), ")"
)
)
}
if (!is.null(x$bias$trimfill)) {
adj <- .r4vn_meta_prediction(
x$bias$trimfill, x$config, x$ci, x$analysis_data$.n_for_back
)
bias_bits <- c(
bias_bits,
paste0(
"trim-and-fill imputed ", x$bias$trimfill$k0 %||% 0L,
" study/studies and gave an adjusted pooled estimate of ",
ci(adj$estimate, adj$lower, adj$upper)
)
)
}
if (length(bias_bits)) {
add(
"Small-study effects / publication bias",
paste0(
paste(bias_bits, collapse = "; "),
". Funnel asymmetry and trim-and-fill are sensitivity diagnostics and do not by themselves prove or exclude publication bias."
)
)
} else if (!is.null(x$bias$status)) {
add("Small-study effects / publication bias", x$bias$status)
}
}
if (!is.null(x$influence$table)) {
influential <- x$influence$table$Study[
x$influence$table$Influential %in% "Yes"
]
add(
"Influence diagnostics",
if (length(influential)) {
paste0(
"The following studies were flagged as influential: ",
paste(influential, collapse = ", "),
". Their data and analytic assumptions should be reviewed."
)
} else {
"No study was flagged as influential by the fitted influence diagnostics."
}
)
}
if (!is.null(x$leaveout$estimate) && length(x$leaveout$estimate)) {
rng <- range(x$leaveout$estimate, finite = TRUE)
add(
"Leave-one-out analysis",
paste0(
"Across leave-one-out analyses, pooled estimates ranged from ",
.r4vn_meta_fmt(rng[1L], x$digits), " to ",
.r4vn_meta_fmt(rng[2L], x$digits),
". Compare this range with the main pooled estimate to assess robustness to individual studies."
)
)
}
if (!is.null(x$cumulative$plot_data) && nrow(x$cumulative$plot_data)) {
z <- x$cumulative$plot_data
add(
"Cumulative meta-analysis",
paste0(
"The cumulative pooled estimate changed from ",
ci(z$Estimate[1L], z$Lower[1L], z$Upper[1L]), " at ",
x$cumulative$label %||% x$cumulative$order, " = ", z[[1L]][1L],
" to ", ci(z$Estimate[nrow(z)], z$Lower[nrow(z)], z$Upper[nrow(z)]),
" at ", x$cumulative$label %||% x$cumulative$order, " = ",
z[[1L]][nrow(z)],
". The cumulative table and plot show how evidence evolved as studies accumulated."
)
)
}
if (!length(rows)) {
return(data.frame(Section = character(), Interpretation = character(),
stringsAsFactors = FALSE))
}
out <- do.call(rbind, rows)
rownames(out) <- NULL
out
}
.r4vn_meta_results_text <- function(x) {
z <- .r4vn_meta_interpretation(x)
paste(z$Interpretation, collapse = " ")
}
.r4vn_meta_html_table <- function(data, title = NULL) {
if (is.null(data) || !is.data.frame(data)) return("")
th <- paste0("<th>", .r4vn_meta_escape(names(data)), "</th>",
collapse = "")
body <- if (!nrow(data)) "" else paste(
vapply(seq_len(nrow(data)), function(i) {
vals <- vapply(data[i, , drop = FALSE], function(z) {
if (is.na(z)) "" else as.character(z)
}, character(1))
paste0(
"<tr>",
paste0("<td>", .r4vn_meta_escape(vals), "</td>", collapse = ""),
"</tr>"
)
}, character(1)),
collapse = ""
)
heading <- if (is.null(title) || !nzchar(title)) "" else
paste0("<h2>", .r4vn_meta_escape(title), "</h2>")
paste0(
heading,
"<div class=\"table-wrap\"><table><thead><tr>", th,
"</tr></thead><tbody>", body, "</tbody></table></div>"
)
}
.r4vn_meta_report_html <- function(x) {
sections <- character()
sections <- c(sections, .r4vn_meta_html_table(x$overview, "Analysis overview"))
sections <- c(sections, .r4vn_meta_html_table(x$table, "Main analysis"))
sections <- c(
sections,
.r4vn_meta_html_table(x$significance_tests, "Statistical significance tests")
)
sections <- c(
sections,
.r4vn_meta_html_table(x$heterogeneity, "Model and heterogeneity")
)
if (length(x$subgroup_tables)) {
for (nm in names(x$subgroup_tables)) {
subgroup_label <- x$subgroups[[nm]]$label %||% nm
sections <- c(
sections,
.r4vn_meta_html_table(
x$subgroup_tables[[nm]],
paste0("Subgroup analysis: ", subgroup_label)
)
)
}
sections <- c(
sections,
.r4vn_meta_html_table(
x$tables$Subgroup_test, "Test for subgroup differences"
)
)
}
if (!is.null(x$regression$table)) {
sections <- c(
sections,
.r4vn_meta_html_table(x$regression$table, "Meta-regression"),
.r4vn_meta_html_table(x$regression$statistics,
"Meta-regression statistics")
)
if (!is.null(x$regression$univariable)) {
sections <- c(
sections,
.r4vn_meta_html_table(x$regression$univariable,
"Univariable moderator screening")
)
}
}
if (!is.null(x$bias$table)) {
sections <- c(
sections,
.r4vn_meta_html_table(
x$bias$table, "Small-study effects / publication bias"
)
)
}
if (!is.null(x$bias$status)) {
sections <- c(
sections,
paste0("<p class=\"note\">", .r4vn_meta_escape(x$bias$status), "</p>")
)
}
if (!is.null(x$bias$adjusted)) {
sections <- c(
sections,
.r4vn_meta_html_table(x$bias$adjusted,
"Publication-bias sensitivity estimates")
)
}
if (!is.null(x$small_sample$comparison)) {
sections <- c(
sections,
.r4vn_meta_html_table(x$small_sample$comparison,
"Few-study inference sensitivity"),
paste0("<p class=\"note\">",
.r4vn_meta_escape(x$small_sample$status), "</p>")
)
}
if (!is.null(x$leaveout$table)) {
sections <- c(
sections,
.r4vn_meta_html_table(
x$leaveout$table, "Leave-one-out sensitivity analysis"
)
)
}
if (!is.null(x$influence$table)) {
sections <- c(
sections,
.r4vn_meta_html_table(x$influence$table, "Influence diagnostics")
)
}
if (!is.null(x$cumulative$table)) {
sections <- c(
sections,
.r4vn_meta_html_table(
x$cumulative$table, "Cumulative meta-analysis"
)
)
}
if (isTRUE(x$report)) {
narrative <- if (!is.null(x$interpretation) && nrow(x$interpretation)) {
paste(
vapply(seq_len(nrow(x$interpretation)), function(i) {
paste0(
"<h3>", .r4vn_meta_escape(x$interpretation$Section[i]), "</h3>",
"<p>", .r4vn_meta_escape(x$interpretation$Interpretation[i]), "</p>"
)
}, character(1)),
collapse = "\n"
)
} else {
paste0("<p>", .r4vn_meta_escape(x$results_text), "</p>")
}
sections <- c(sections, paste0("<h2>Interpretation</h2>", narrative))
}
if (length(x$plots)) {
for (nm in names(x$plots)) {
sections <- c(
sections,
paste0(
"<h2>", .r4vn_meta_escape(x$plot_titles[[nm]]), "</h2>",
"<div class=\"figure\"><img src=\"",
.r4vn_meta_escape(basename(x$plots[[nm]])),
"\" alt=\"", .r4vn_meta_escape(x$plot_titles[[nm]]),
"\"></div>"
)
)
}
}
css <- paste0(
"body{font-family:'Times New Roman',Times,serif;color:#111;background:#fff;",
"margin:24px;line-height:1.35}",
".report{max-width:1180px;margin:auto}",
"h1{font-size:22px;margin:0 0 16px}",
"h2{font-size:18px;margin:24px 0 8px}",
"h3{font-size:15px;margin:15px 0 5px}",
".table-wrap{overflow-x:auto;margin-bottom:14px}",
"table{border-collapse:collapse;width:100%;border-top:2px solid #111;",
"border-bottom:2px solid #111}",
"th{border-bottom:1.5px solid #111;padding:6px 8px;text-align:center;",
"font-weight:700;white-space:nowrap}",
"td{padding:5px 8px;text-align:center;vertical-align:top}",
"th:first-child,td:first-child{text-align:left}",
".note{font-size:13px}",
".figure{text-align:center;margin:8px 0 22px}",
".figure img{max-width:100%;height:auto}",
".foot{font-size:12px;margin-top:18px;color:#444}"
)
paste0(
"<!doctype html><html><head><meta charset=\"utf-8\">",
"<meta name=\"viewport\" content=\"width=device-width,initial-scale=1\">",
"<title>", .r4vn_meta_escape(x$title), "</title><style>", css,
"</style></head><body><main class=\"report\"><h1>",
.r4vn_meta_escape(x$title), "</h1>",
paste(sections, collapse = "\n"),
"<div class=\"foot\">R4VN meta-analysis. Random-effects estimator: ",
.r4vn_meta_escape(x$method), "; inference: ",
.r4vn_meta_escape(x$test_method), ".",
"</div></main></body></html>"
)
}
.r4vn_meta_open <- function(path) {
viewer <- getOption("viewer")
p <- normalizePath(path, winslash = "/", mustWork = TRUE)
if (is.function(viewer)) viewer(p) else utils::browseURL(p)
invisible(path)
}
#' Publication-Ready Meta-Analysis in One Command
#'
#' Fits fixed- and/or random-effects meta-analysis and returns a complete,
#' publication-ready report. `tabmeta()` accepts a binary 2-by-2 table
#' (`a`, `b`, `c`, `d`), event/total data, continuous summaries, rates,
#' correlations, or study-level effect estimates with standard errors or
#' confidence intervals. The default `profile="auto"` adds prediction,
#' few-study inference, subgroup/moderator results when requested,
#' small-study-effect diagnostics, sensitivity analyses, and figures.
#'
#' @param data Optional data frame. When omitted, active R4VN data are used.
#' @param study Study label variable.
#' @param effect Generic study-level effect estimate on the natural scale.
#' @param se Standard error on the analysis scale. For ratio measures this is
#' the standard error of the log effect.
#' @param lower,upper Lower and upper confidence limits for `effect`.
#' @param a,b,c,d Binary 2-by-2 cells: events and non-events in group 1
#' (`a`, `b`) and group 0 (`c`, `d`). Supplying these four arguments is
#' equivalent to `event1=a`, `n1=a+b`, `event0=c`, `n0=c+d`.
#' @param event1,n1,event0,n0 Events and total sample sizes in groups 1 and 0.
#' @param mean1,sd1,mean0,sd0 Group means and standard deviations.
#' @param event,n,time Single-group events, sample size, and person-time.
#' @param time1,time0 Person-time in groups 1 and 0 for incidence-rate ratios.
#' @param or,rr,rd,hr,irr,md,smd,prop,rate,cor Logical effect selectors.
#' With raw binary input, OR is inferred if none is selected; with continuous
#' summaries, MD is inferred. Set a selector to request another measure.
#' @param fixed Fit a fixed-effect model in addition to, or instead of, the
#' random-effects model.
#' @param random Fit a random-effects model; default `TRUE`.
#' @param method Random-effects tau-squared estimator; default `"REML"`.
#' @param hk Backward-compatible logical shortcut. `TRUE` is
#' `small="knha"`; `FALSE` is `small="z"`. Prefer `small` in new code.
#' @param small Random-effects inference: `"auto"`, `"adhoc"`, `"knha"`,
#' `"t"`, or `"z"`. `"auto"` uses modified Hartung-Knapp (`"adhoc"`)
#' when there are at most 10 studies and the normal approximation otherwise.
#' @param prediction Add a prediction interval. `NULL` uses the profile default.
#' @param by,subgroup One or more subgroup variables. Use one unquoted variable,
#' a character vector, or `vars(region, design)`. `subgroup` is a readable
#' alias for `by`; both may be combined. Factor levels and value labels from
#' labelled data are used in tables, interpretation, and subgroup figures.
#' For every estimable subgroup, the result includes its pooled effect and
#' confidence interval, effect test/df/p-value, prediction interval,
#' Cochran's Q/df/p-value, I-squared and tau-squared with confidence
#' intervals when estimable. The publication table is transposed: statistics
#' are rows and subgroup labels are columns. Multiple subgroup variables
#' produce one transposed table per variable.
#' @param reg,moderator One or more meta-regression moderators created with
#' `vars()` or supplied as character names. `moderator` is an alias for
#' `reg`. The output includes the multivariable model and a univariable
#' moderator screen.
#' @param bias Assess small-study effects/publication bias. `NULL` uses the
#' profile default.
#' @param bias_methods One or more of `"egger"`, `"begg"`, `"trimfill"`,
#' `"failsafe"`, or `"selection"`. `"auto"` runs Egger, Begg, and
#' trim-and-fill; the `"full"` profile also requests fail-safe N and a
#' selection model. Methods are used as sensitivity diagnostics, not as proof
#' that publication bias is or is not present.
#' @param leaveout Perform leave-one-out sensitivity analysis. `NULL` uses the
#' profile default.
#' @param influence Perform influence diagnostics. `NULL` uses the profile
#' default.
#' @param cumulative Optional variable defining the ordering for cumulative
#' meta-analysis, usually publication year.
#' @param transform Single-proportion transformation: `"logit"`, `"arcsine"`,
#' `"ft"`, or `"none"`.
#' @param cc Continuity correction for zero events; default `0.5`.
#' @param zero Handling of double-zero binary studies: `"keep"` or `"exclude"`.
#' @param ci Confidence level as a proportion.
#' @param digit Number of decimals for effect estimates and confidence limits,
#' including subgroup and meta-regression estimates. Default `2`.
#' @param p_digit Number of decimals for p-values.
#' @param profile Analysis profile. `"auto"` (default) creates a complete
#' appropriate analysis; `"brief"` creates the main model and forest plot;
#' `"full"` also requests extended bias diagnostics and radial/L'Abbe plots;
#' `"custom"` turns optional modules off unless explicitly requested.
#' @param full Backward-compatible profile shortcut: `TRUE` selects `"full"`
#' and `FALSE` selects `"custom"`.
#' @param plot Create publication-ready plots. `NULL` uses the profile default.
#' With subgroup analysis, the HTML Viewer includes the overall forest plot
#' and one clearly titled forest plot for every estimable subgroup level.
#' @param plot_display Plot type(s) also drawn in the interactive R/RStudio Plot
#' pane when `plot=TRUE`. Default `"forest"`. Use a character vector for Plot
#' history, `"all"` for every available figure, or `"none"` to suppress Plot
#' pane drawing while retaining figures in the HTML Viewer. When subgroup
#' analysis is present, selecting `"forest"` adds the overall forest plot and
#' every subgroup forest plot to Plot history; use Previous/Next to review.
#' @param plot_args Named list of plot options. Supply common options directly,
#' or nested lists such as `list(all=list(color="navy"),
#' forest=list(xlim=c(0.2, 2)))`.
#' @param report Add manuscript-style results text.
#' @param interpretation Add a detailed, sectioned interpretation covering the
#' pooled effect, heterogeneity, prediction interval, subgroup differences,
#' moderators, small-study effects, few-study inference, influence,
#' leave-one-out robustness, and cumulative evidence whenever available.
#' Default is `FALSE`.
#' @param title Report title.
#' @param show Open the complete HTML report; default `TRUE`.
#' @param export Optional direct export format(s): `"html"`, `"docx"`,
#' `"xlsx"`, `"pdf"`, or `"png"`.
#' @param file Export file or base path. The extension may determine the format.
#' @param open Open the exported file.
#' @param strict Stop when an optional diagnostic or figure cannot be created.
#' The default `FALSE` retains the main analysis and issues a warning/note.
#'
#' @return Invisibly returns an object of class `r4vn_meta` and `r4vn_tab`.
#' Stable publication components are available in `$estimates`, `$tests`,
#' `$diagnostics`, `$models`, `$tables`, and `$metadata`. The table
#' `$tables$Statistical_tests` gives the pooled-effect and Cochran's Q tests.
#' Its significance and conclusion columns are added only when
#' `interpretation=TRUE`. `$tables$Subgroup` is the transposed publication
#' table for the first subgroup variable, `$subgroup_tables` contains every
#' transposed subgroup table, and `$subgroup_long` retains tidy long output.
#' @family R4VN tables
#' @seealso `tabexport`, `vars`
#' @export
#'
#' @examples
#' if (requireNamespace("metafor", quietly = TRUE)) {
#' dat <- read.csv(
#' system.file("extdata", "meta_example.csv", package = "R4VN")
#' )
#'
#' # 1. Binary outcome from a, b, c, d. OR is inferred automatically.
#' dat$non_event_treat <- dat$n_treat - dat$event_treat
#' dat$non_event_control <- dat$n_control - dat$event_control
#' m_abcd <- tabmeta(
#' data = dat, study = study,
#' a = event_treat, b = non_event_treat,
#' c = event_control, d = non_event_control,
#' profile = "custom", plot = FALSE, show = FALSE
#' )
#' m_abcd$estimates$overall
#' m_abcd$tables$Binary_2x2
#' m_abcd$tables$Statistical_tests
#' m_abcd$tests$overall
#' m_abcd$tests$heterogeneity
#' plot(
#' m_abcd, type = "forest", show_abcd = TRUE,
#' abcd_titles = c("Events T", "No event T", "Events C", "No event C")
#' )
#'
#' # 2. Equivalent event/total syntax; request RR or RD with rr/rd=TRUE.
#' m_or <- tabmeta(
#' data = dat, study = study,
#' event1 = event_treat, n1 = n_treat,
#' event0 = event_control, n0 = n_control,
#' or = TRUE, profile = "custom", plot = FALSE, show = FALSE
#' )
#'
#' # 3. Complete automatic analysis: report, prediction, diagnostics, plots.
#' \donttest{
#' m_all <- tabmeta(
#' data = dat, study = study,
#' a = event_treat, b = non_event_treat,
#' c = event_control, d = non_event_control,
#' show = FALSE
#' )
#' }
#'
#' # 4. Subgroup analysis: one or several subgroup variables.
#' m_sub <- tabmeta(
#' data = dat, study = study,
#' effect = OR, lower = LCI, upper = UCI, or = TRUE,
#' subgroup = region,
#' profile = "custom", plot = FALSE, show = FALSE
#' )
#' m_sub$tables$Subgroup
#' m_sub$subgroup_long
#' m_sub$tables$Subgroup_test
#' # Each subgroup row includes its effect test and heterogeneity test.
#' names(m_sub$tables$Subgroup)
#' # If the subgroup variable carries value labels, R4VN prints those labels
#' # instead of numeric codes. With plot=TRUE, both Viewer and Plot history
#' # contain the overall plot and one clearly titled plot per subgroup level.
#' dat$risk_group <- rep(c(1, 2), length.out = nrow(dat))
#' attr(dat$risk_group, "label") <- "Baseline risk"
#' attr(dat$risk_group, "labels") <- c("Lower risk" = 1, "Higher risk" = 2)
#' \donttest{
#' m_labelled <- tabmeta(
#' data = dat, study = study,
#' effect = OR, lower = LCI, upper = UCI, or = TRUE,
#' subgroup = risk_group, profile = "custom",
#' plot = TRUE, plot_display = "forest", show = FALSE
#' )
#' m_labelled$plot_titles
#' }
#'
#' # Multiple subgroup analyses can be requested together.
#' dat$period <- ifelse(dat$year < median(dat$year), "Earlier", "Later")
#' m_sub2 <- tabmeta(
#' data = dat, study = study,
#' effect = OR, lower = LCI, upper = UCI, or = TRUE,
#' subgroup = vars(region, period),
#' profile = "custom", plot = FALSE, show = FALSE
#' )
#'
#' # 5. Meta-regression with numeric and categorical moderators.
#' attr(dat$year, "label") <- "Publication year"
#' attr(dat$region, "label") <- "Geographic region"
#' m_reg <- tabmeta(
#' data = dat, study = study,
#' effect = OR, lower = LCI, upper = UCI, or = TRUE,
#' moderator = vars(c.year, region),
#' profile = "custom", plot = FALSE, show = FALSE
#' )
#' m_reg$tables$Meta_regression
#' m_reg$tables$Moderator_univariable
#' # Intercept is written in full; moderator labels are used when available.
#' # Set digit=4, for example, when four decimal places are required.
#'
#' # 6. Publication-bias and small-study-effect sensitivity analyses.
#' m_bias <- tabmeta(
#' data = dat, study = study,
#' effect = OR, lower = LCI, upper = UCI, or = TRUE,
#' bias = TRUE,
#' bias_methods = c("egger", "begg", "trimfill", "failsafe"),
#' profile = "custom", plot = FALSE, show = FALSE
#' )
#' m_bias$tables$Publication_bias
#' m_bias$tables$Publication_bias_adjusted
#'
#' # 7. Few-study inference. auto uses modified Hartung-Knapp at <=10 studies.
#' m_few <- tabmeta(
#' data = dat[1:8, ], study = study,
#' effect = OR, lower = LCI, upper = UCI, or = TRUE,
#' small = "auto", profile = "custom", plot = FALSE, show = FALSE
#' )
#' m_few$tables$Small_sample_inference
#'
#' # 8. Generic hazard ratios with confidence intervals.
#' m_hr <- tabmeta(
#' data = dat, study = study,
#' effect = OR, lower = LCI, upper = UCI, hr = TRUE,
#' profile = "custom", plot = FALSE, show = FALSE
#' )
#'
#' # 9. Continuous MD/SMD, single proportion/rate, correlation, and IRR.
#' cont <- data.frame(
#' study = paste0("C", 1:5),
#' m1 = c(12, 14, 13, 16, 15), s1 = c(3, 4, 3, 5, 4), n1 = rep(60, 5),
#' m0 = c(15, 15, 16, 18, 16), s0 = c(4, 4, 5, 5, 4), n0 = rep(60, 5)
#' )
#' m_md <- tabmeta(
#' cont, study, mean1 = m1, sd1 = s1, n1 = n1,
#' mean0 = m0, sd0 = s0, n0 = n0, md = TRUE,
#' profile = "custom", plot = FALSE, show = FALSE
#' )
#' m_smd <- tabmeta(
#' cont, study, mean1 = m1, sd1 = s1, n1 = n1,
#' mean0 = m0, sd0 = s0, n0 = n0, smd = TRUE,
#' profile = "custom", plot = FALSE, show = FALSE
#' )
#'
#' one <- data.frame(
#' study = paste0("P", 1:5), events = c(8, 12, 15, 10, 14),
#' total = c(100, 110, 120, 90, 105), person_time = c(80, 90, 95, 75, 88),
#' correlation = c(.20, .28, .15, .31, .24)
#' )
#' m_prop <- tabmeta(
#' one, study, event = events, n = total, prop = TRUE,
#' profile = "custom", plot = FALSE, show = FALSE
#' )
#' m_rate <- tabmeta(
#' one, study, event = events, time = person_time, rate = TRUE,
#' profile = "custom", plot = FALSE, show = FALSE
#' )
#' m_cor <- tabmeta(
#' one, study, effect = correlation, n = total, cor = TRUE,
#' profile = "custom", plot = FALSE, show = FALSE
#' )
#' m_irr <- tabmeta(
#' dat[1:5, ], study,
#' event1 = event_treat, time1 = n_treat,
#' event0 = event_control, time0 = n_control, irr = TRUE,
#' profile = "custom", plot = FALSE, show = FALSE
#' )
#'
#' # 10. Cumulative meta-analysis ordered by publication year.
#' m_cum <- tabmeta(
#' data = dat, study = study,
#' effect = OR, lower = LCI, upper = UCI, or = TRUE,
#' cumulative = year,
#' profile = "custom", plot = FALSE, show = FALSE
#' )
#' m_cum$tables$Cumulative
#'
#' # 11. Draw or save individual publication figures.
#' if (interactive()) {
#' # plot=TRUE keeps all figures in the HTML Viewer and also draws the
#' # selected figures in the R/RStudio Plot pane and Plot history.
#' m_publication <- tabmeta(
#' data = dat, study = study,
#' a = event_treat, b = non_event_treat,
#' c = event_control, d = non_event_control,
#' interpretation = TRUE,
#' plot = TRUE,
#' plot_display = c("forest", "funnel", "trimfill"),
#' plot_args = list(
#' all = list(
#' font_family = "Arial", background = "white",
#' title_color = "#17365D", title_size = 1.1,
#' text_size = 0.86, axis_size = 0.92
#' ),
#' forest = list(
#' subtitle = "Random-effects model with 95% confidence intervals",
#' caption = "Square size reflects study weight; diamond is pooled effect.",
#' margins = c(5.5, 4.2, 5.0, 2.0),
#' point_color = "#1F4E79", ci_color = "#5B9BD5",
#' summary_color = "#C00000", summary_border = "#7F0000",
#' point_shape = 15, row_shade = "zebra",
#' shade_color = "#F5F7FA", show_weights = TRUE,
#' weight_title = "Weight", estimate_title = "OR (95% CI)",
#' show_prediction = FALSE,
#' ref_color = "#666666", ref_type = 2,
#' xlim = c(0.2, 2.0), ticks = c(0.25, 0.5, 1, 1.5, 2)
#' ),
#' funnel = list(
#' point_shape = 21, point_color = "#1F4E79",
#' point_bg = "#D9EAF7", point_size = 1.1,
#' contour_levels = c(90, 95, 99),
#' contour_colors = c("#FFF2CC", "#FCE4D6", "#E2F0D9"),
#' funnel_label = "out", funnel_legend = "topright"
#' )
#' )
#' )
#' m_publication$tables$Interpretation
#'
#' # Any figure can be redrawn or saved independently.
#' plot(m_publication, type = "forest")
#' plot(m_publication, type = "funnel", contour = TRUE)
#' plot(m_publication, type = "trimfill", contour = TRUE)
#' plot(m_reg, type = "bubble", moderator = "year")
#' plot(
#' m_publication, type = "forest", file = "forest_publication.tiff",
#' width = 2400, height = 1800, res = 300,
#' font_family = "Arial", point_color = "#1F4E79",
#' ci_color = "#5B9BD5", summary_color = "#C00000",
#' show_weights = TRUE, show_prediction = FALSE,
#' xlim = c(0.5, 1.5), ticks = c(0.5, 0.75, 1, 1.25, 1.5)
#' )
#' # Values outside xlim remain exact in the Estimate (95% CI) column;
#' # the graphical confidence interval is clipped with an arrow.
#' m_publication$overall[c("pi_lower", "pi_upper")]
#' tabmeta(
#' data = dat, study = study,
#' a = event_treat, b = non_event_treat,
#' c = event_control, d = non_event_control,
#' export = c("docx", "xlsx"), file = "meta_report",
#' show = FALSE
#' )
#' }
#' }
tabmeta <- function(
data = NULL,
study,
effect = NULL, se = NULL, lower = NULL, upper = NULL,
a = NULL, b = NULL, c = NULL, d = NULL,
event1 = NULL, n1 = NULL, event0 = NULL, n0 = NULL,
mean1 = NULL, sd1 = NULL, mean0 = NULL, sd0 = NULL,
event = NULL, n = NULL, time = NULL, time1 = NULL, time0 = NULL,
or = FALSE, rr = FALSE, rd = FALSE, hr = FALSE, irr = FALSE,
md = FALSE, smd = FALSE, prop = FALSE, rate = FALSE, cor = FALSE,
fixed = FALSE, random = TRUE, method = "REML", hk = NULL,
small = c("auto", "adhoc", "knha", "t", "z"),
prediction = NULL, by = NULL, subgroup = NULL,
reg = NULL, moderator = NULL,
bias = NULL, bias_methods = "auto",
leaveout = NULL, influence = NULL, cumulative = NULL,
transform = NULL, cc = 0.5, zero = c("keep", "exclude"),
ci = 0.95, digit = 2, p_digit = 3,
profile = c("auto", "brief", "full", "custom"),
full = NULL, plot = NULL, plot_display = "forest", plot_args = list(),
report = FALSE, interpretation = FALSE,
title = NULL, show = TRUE,
export = NULL, file = NULL, open = FALSE, strict = FALSE) {
# `c` is a public 2-by-2 cell argument. Capture its expression before
# restoring base::c locally; otherwise the argument promise can mask c()
# and be forced while ordinary vectors are being constructed.
c_expr <- substitute(c)
c <- base::c
.r4vn_meta_need()
data <- .r4vn_resolve_analysis_data(data)
env <- parent.frame()
settings <- .r4vn_meta_profile(
profile = match.arg(profile), full = full,
prediction = prediction, bias = bias, leaveout = leaveout,
influence = influence, plot = plot
)
profile <- settings$profile
prediction <- settings$prediction
bias <- settings$bias
leaveout <- settings$leaveout
influence <- settings$influence
plot <- settings$plot
full <- settings$full
flag_names <- c(
"or", "rr", "rd", "hr", "irr", "md", "smd", "prop", "rate", "cor",
"fixed", "random", "prediction", "bias", "leaveout",
"influence", "full", "plot", "report", "interpretation", "show",
"open", "strict"
)
for (z in flag_names) .r4vn_meta_flag(get(z), z)
if (!is.null(hk)) .r4vn_meta_flag(hk, "hk")
if (!is.list(plot_args) ||
(length(plot_args) && is.null(names(plot_args)))) {
stop("`plot_args` must be a named list.", call. = FALSE)
}
plot_types <- c(
"forest", "funnel", "trimfill", "leaveout", "influence", "baujat",
"radial", "labbe", "cumulative", "bubble"
)
if (is.null(plot_display) || identical(plot_display, FALSE)) {
plot_display <- character()
} else {
plot_display <- unique(tolower(as.character(plot_display)))
plot_display <- plot_display[!is.na(plot_display) & nzchar(plot_display)]
if ("none" %in% plot_display) plot_display <- character()
bad_display <- setdiff(plot_display, c(plot_types, "all"))
if (length(bad_display)) {
stop("Unsupported `plot_display`: ", paste(bad_display, collapse = ", "),
".", call. = FALSE)
}
}
if (!fixed && !random) {
stop("At least one of `fixed` or `random` must be TRUE.", call. = FALSE)
}
.r4vn_meta_num1(ci, "ci", 0, 1, FALSE, FALSE)
.r4vn_meta_num1(cc, "cc", 0, Inf)
digit <- as.integer(.r4vn_meta_num1(digit, "digit", 0, 10))
p_digit <- as.integer(.r4vn_meta_num1(p_digit, "p_digit", 0, 10))
method_input <- toupper(as.character(method)[1L])
method_map <- c(
REML = "REML", ML = "ML", DL = "DL", HE = "HE", HS = "HS",
HSK = "HSk", SJ = "SJ", EB = "EB", PM = "PM", PMM = "PMM",
GENQ = "GENQ"
)
if (!method_input %in% names(method_map)) {
stop("Unsupported random-effects method: ", method_input, ".",
call. = FALSE)
}
method <- unname(method_map[[method_input]])
zero <- match.arg(zero)
small <- match.arg(tolower(as.character(small)[1L]),
c("auto", "adhoc", "knha", "t", "z"))
if (!is.null(hk)) small <- if (isTRUE(hk)) "knha" else "z"
bias_methods <- .r4vn_meta_bias_methods(bias_methods, profile)
report_text <- isTRUE(report) || isTRUE(interpretation)
expr <- list(
study = substitute(study),
effect = substitute(effect), se = substitute(se),
lower = substitute(lower), upper = substitute(upper),
a = substitute(a), b = substitute(b), c = c_expr, d = substitute(d),
event1 = substitute(event1), n1 = substitute(n1),
event0 = substitute(event0), n0 = substitute(n0),
mean1 = substitute(mean1), sd1 = substitute(sd1),
mean0 = substitute(mean0), sd0 = substitute(sd0),
event = substitute(event), n = substitute(n),
time = substitute(time), time1 = substitute(time1), time0 = substitute(time0),
by = substitute(by), subgroup = substitute(subgroup),
cumulative = substitute(cumulative)
)
study_value <- .r4vn_meta_eval(expr$study, data, env, "study", TRUE)
if (length(study_value) != nrow(data)) {
stop("`study` must have one value per row of `data`.", call. = FALSE)
}
study_value <- as.character(study_value)
if (anyNA(study_value) || any(!nzchar(study_value))) {
stop("`study` cannot contain missing or empty labels.", call. = FALSE)
}
study_value <- make.unique(study_value)
effect_value <- .r4vn_meta_eval(expr$effect, data, env, "effect")
getv <- function(nm) .r4vn_meta_eval(expr[[nm]], data, env, nm)
values <- lapply(
c("se", "lower", "upper", "a", "b", "c", "d",
"event1", "n1", "event0", "n0",
"mean1", "sd1", "mean0", "sd0", "event", "n", "time",
"time1", "time0"),
getv
)
names(values) <- c(
"se", "lower", "upper", "a", "b", "c", "d",
"event1", "n1", "event0", "n0",
"mean1", "sd1", "mean0", "sd0", "event", "n", "time",
"time1", "time0"
)
values$effect <- effect_value
for (nm in names(values)) {
if (!is.null(values[[nm]]) && length(values[[nm]]) != nrow(data)) {
stop("`", nm, "` must have one value per row of `data`.",
call. = FALSE)
}
}
abcd_present <- vapply(values[c("a", "b", "c", "d")],
Negate(is.null), logical(1))
if (any(abcd_present) && !all(abcd_present)) {
stop("Supply all four binary cells: `a`, `b`, `c`, and `d`.",
call. = FALSE)
}
if (all(abcd_present)) {
old_binary <- vapply(values[c("event1", "n1", "event0", "n0")],
Negate(is.null), logical(1))
if (any(old_binary)) {
stop("Use either `a,b,c,d` or `event1,n1,event0,n0`, not both.",
call. = FALSE)
}
values$event1 <- as.numeric(values$a)
values$n1 <- as.numeric(values$a) + as.numeric(values$b)
values$event0 <- as.numeric(values$c)
values$n0 <- as.numeric(values$c) + as.numeric(values$d)
}
inferred <- if (all(abcd_present) ||
all(vapply(values[c("event1", "n1", "event0", "n0")],
Negate(is.null), logical(1)))) {
"OR"
} else if (all(vapply(values[c("event1", "time1", "event0", "time0")],
Negate(is.null), logical(1)))) {
"IRR"
} else if (all(vapply(values[c("mean1", "sd1", "n1", "mean0", "sd0", "n0")],
Negate(is.null), logical(1)))) {
"MD"
} else if (all(vapply(values[c("event", "n")],
Negate(is.null), logical(1)))) {
"PROP"
} else if (all(vapply(values[c("event", "time")],
Negate(is.null), logical(1)))) {
"RATE"
} else NULL
type <- .r4vn_meta_effect_type(
or, rr, rd, hr, irr, md, smd, prop, rate, cor,
!is.null(effect_value), inferred = inferred
)
config <- .r4vn_meta_config(type, transform)
if (identical(config$scale, "ft")) {
warning(
paste(
"Freeman-Tukey transformed proportions are supported, but",
"`transform=\"logit\"` is the R4VN default and is generally preferred."
),
call. = FALSE
)
}
by_names <- unique(c(
.r4vn_meta_varnames(expr$by, data, env, "by", TRUE),
.r4vn_meta_varnames(expr$subgroup, data, env, "subgroup", TRUE)
))
by_name <- if (length(by_names)) by_names[1L] else NULL
cumulative_name <- .r4vn_meta_varname(
expr$cumulative, data, env, "cumulative", TRUE
)
reg_names <- unique(c(.r4vn_meta_reg_names(reg, data),
.r4vn_meta_reg_names(moderator, data)))
raw_required <- switch(
type,
OR = c("event1", "n1", "event0", "n0"),
RR = c("event1", "n1", "event0", "n0"),
RD = c("event1", "n1", "event0", "n0"),
IRR = c("event1", "time1", "event0", "time0"),
MD = c("mean1", "sd1", "n1", "mean0", "sd0", "n0"),
SMD = c("mean1", "sd1", "n1", "mean0", "sd0", "n0"),
PROP = c("event", "n"),
RATE = c("event", "time"),
COR = c("effect", "n"),
HR = character(),
GENERIC = character()
)
use_generic <- !is.null(effect_value) &&
(type %in% c("HR", "OR", "RR", "IRR", "RD", "MD", "SMD", "GENERIC") ||
(type == "COR" && is.null(values$n)) ||
type %in% c("PROP", "RATE"))
if (!use_generic) {
absent <- raw_required[
vapply(raw_required, function(z) is.null(values[[z]]), logical(1))
]
if (length(absent)) {
stop("Missing required input: ", paste(absent, collapse = ", "), ".",
call. = FALSE)
}
}
if (use_generic) {
yi <- .r4vn_meta_transform(as.numeric(effect_value), config)
if (!is.null(values$se)) {
sei <- as.numeric(values$se)
if (any(sei <= 0, na.rm = TRUE)) {
stop("`se` must be greater than 0.", call. = FALSE)
}
} else if (!is.null(values$lower) && !is.null(values$upper)) {
lower_t <- .r4vn_meta_transform(as.numeric(values$lower), config)
upper_t <- .r4vn_meta_transform(as.numeric(values$upper), config)
if (any(lower_t >= upper_t, na.rm = TRUE)) {
stop("Each `lower` value must be below its `upper` value.",
call. = FALSE)
}
zcrit <- stats::qnorm(1 - (1 - ci) / 2)
sei <- (upper_t - lower_t) / (2 * zcrit)
} else {
stop("Generic `effect` input requires `se` or both `lower` and `upper`.",
call. = FALSE)
}
vi <- sei^2
} else {
es <- .r4vn_meta_make_es(type, config, values, cc, zero)
if (is.null(es)) {
stop("This input combination is not supported.", call. = FALSE)
}
yi <- as.numeric(es$yi)
vi <- as.numeric(es$vi)
}
analysis <- data.frame(
.row = seq_len(nrow(data)),
.study = study_value,
yi = yi,
vi = vi,
stringsAsFactors = FALSE
)
if (type %in% c("PROP", "COR") && !is.null(values$n)) {
analysis$.n_for_back <- as.numeric(values$n)
} else {
analysis$.n_for_back <- rep(NA_real_, nrow(data))
}
extra_names <- unique(c(by_names, cumulative_name, reg_names))
extra_names <- extra_names[!is.na(extra_names) & nzchar(extra_names)]
if (length(extra_names)) {
for (nm in extra_names) analysis[[nm]] <- data[[nm]]
}
if (!is.null(values$event1)) analysis$.event1 <- as.numeric(values$event1)
if (!is.null(values$n1)) analysis$.n1 <- as.numeric(values$n1)
if (!is.null(values$event0)) analysis$.event0 <- as.numeric(values$event0)
if (!is.null(values$n0)) analysis$.n0 <- as.numeric(values$n0)
if (!is.null(values$mean1)) analysis$.mean1 <- as.numeric(values$mean1)
if (!is.null(values$sd1)) analysis$.sd1 <- as.numeric(values$sd1)
if (!is.null(values$mean0)) analysis$.mean0 <- as.numeric(values$mean0)
if (!is.null(values$sd0)) analysis$.sd0 <- as.numeric(values$sd0)
if (!is.null(values$event)) analysis$.event <- as.numeric(values$event)
if (!is.null(values$n)) analysis$.n <- as.numeric(values$n)
if (!is.null(values$time)) analysis$.time <- as.numeric(values$time)
if (!is.null(values$time1)) analysis$.time1 <- as.numeric(values$time1)
if (!is.null(values$time0)) analysis$.time0 <- as.numeric(values$time0)
if (all(abcd_present)) {
analysis$.a <- as.numeric(values$a)
analysis$.b <- as.numeric(values$b)
analysis$.c <- as.numeric(values$c)
analysis$.d <- as.numeric(values$d)
}
keep <- is.finite(analysis$yi) & is.finite(analysis$vi) & analysis$vi > 0
if (sum(keep) < 2L) {
stop(
"At least two studies with finite effect estimates and variances are required.",
call. = FALSE
)
}
excluded <- analysis[!keep, c(".study", "yi", "vi"), drop = FALSE]
analysis <- analysis[keep, , drop = FALSE]
rownames(analysis) <- NULL
# Base vector subsetting can drop custom variable/value-label attributes,
# especially when users attach `label` and `labels` without a labelled
# vector class. Rebuild analysis variables from the source data and restore
# these attributes after study exclusion so subgroup output keeps labels.
if (length(extra_names)) {
for (nm in extra_names) {
source_variable <- data[[nm]]
filtered_variable <- source_variable[keep]
for (attribute_name in c("label", "labels", "value.labels")) {
attribute_value <- attr(source_variable, attribute_name, exact = TRUE)
if (!is.null(attribute_value)) {
attr(filtered_variable, attribute_name) <- attribute_value
}
}
analysis[[nm]] <- filtered_variable
}
}
variable_labels <- if (length(extra_names)) {
stats::setNames(
vapply(
extra_names,
function(nm) .r4vn_meta_variable_label(analysis[[nm]], nm),
character(1)
),
extra_names
)
} else character()
test_method <- .r4vn_meta_test_method(small, nrow(analysis), random = random)
random_model <- if (random) {
.r4vn_meta_model(analysis$yi, analysis$vi, method, ci = ci,
fixed = FALSE, test = small)
} else NULL
fixed_model <- if (fixed) {
.r4vn_meta_model(analysis$yi, analysis$vi, method, ci = ci,
fixed = TRUE, test = "z")
} else NULL
primary <- if (!is.null(random_model)) random_model else fixed_model
pred <- .r4vn_meta_prediction(primary, config, ci, analysis$.n_for_back)
overall <- list(
estimate = pred$estimate,
lower = pred$lower,
upper = pred$upper,
pi_lower = pred$pi_lower,
pi_upper = pred$pi_upper,
p = as.numeric(primary$pval)[1L],
statistic = if (!is.null(primary$zval)) {
as.numeric(primary$zval)[1L]
} else NA_real_,
df = if (test_method %in% c("t", "knha", "adhoc")) {
as.numeric((primary$ddf %||% (primary$k - primary$p))[1L])
} else NA_real_,
significant = is.finite(as.numeric(primary$pval)[1L]) &&
as.numeric(primary$pval)[1L] < (1 - ci)
)
subgroup_all <- .r4vn_meta_subgroups_all(
analysis, by_names, method, small, ci, config, digit, p_digit,
fixed = !random, include_interpretation = isTRUE(interpretation)
)
subgroup <- subgroup_all$primary
regression <- .r4vn_meta_regression(
analysis, reg_names, method, small, ci, digit, p_digit, config,
fixed = !random
)
do_prediction <- isTRUE(prediction)
do_bias <- isTRUE(bias)
do_leaveout <- isTRUE(leaveout)
do_influence <- isTRUE(influence)
bias_result <- if (do_bias) {
.r4vn_meta_bias(
primary, primary$k, p_digit, bias_methods, config, ci,
analysis$.n_for_back, strict = strict
)
} else NULL
leave_result <- if (do_leaveout) {
.r4vn_meta_leaveout(primary, config, digit, analysis$.n_for_back)
} else NULL
influence_result <- if (do_influence) {
.r4vn_meta_influence(primary, analysis$.study)
} else NULL
cumulative_result <- .r4vn_meta_cumulative(
analysis, cumulative_name, method, small, ci, config, digit,
fixed = !random
)
small_sample <- if (isTRUE(random)) {
.r4vn_meta_small_sample(
analysis, method, ci, config, analysis$.n_for_back,
requested = small, selected = test_method, digits = digit
)
} else NULL
primary_subgroup <- if (!is.null(by_name)) {
.r4vn_meta_group_values(analysis[[by_name]])
} else NULL
primary_subgroup_title <- if (!is.null(by_name)) {
.r4vn_meta_variable_label(analysis[[by_name]], by_name)
} else "Subgroup"
main_table <- .r4vn_meta_study_table(
analysis, primary, config, digit,
subgroup = primary_subgroup,
subgroup_title = primary_subgroup_title,
ci = ci
)
pooled_rows <- list()
add_pool <- function(label, model) {
pp <- .r4vn_meta_prediction(model, config, ci, analysis$.n_for_back)
row <- as.list(rep("", ncol(main_table)))
names(row) <- names(main_table)
if ("Study" %in% names(row)) row$Study <- label
row[["Effect (95% CI)"]] <- .r4vn_meta_ci_text(
pp$estimate, pp$lower, pp$upper, digit
)
row$Weight <- "100.0%"
as.data.frame(row, stringsAsFactors = FALSE, check.names = FALSE)
}
if (fixed) {
pooled_rows[[length(pooled_rows) + 1L]] <-
add_pool("Overall (fixed)", fixed_model)
}
if (random) {
pooled_rows[[length(pooled_rows) + 1L]] <-
add_pool("Overall (random)", random_model)
}
pool <- do.call(rbind, pooled_rows)
main_table <- rbind(main_table, pool)
rownames(main_table) <- NULL
heterogeneity <- .r4vn_meta_heterogeneity(primary, p_digit)
significance_tests <- .r4vn_meta_significance_tests(
primary, overall, test_method, ci, p_digit,
include_interpretation = isTRUE(interpretation)
)
if (do_prediction &&
is.finite(overall$pi_lower) && is.finite(overall$pi_upper)) {
heterogeneity <- rbind(
heterogeneity,
data.frame(
Statistic = paste0(round(ci * 100), "% prediction interval"),
Value = paste0(
.r4vn_meta_fmt(overall$pi_lower, digit), "-",
.r4vn_meta_fmt(overall$pi_upper, digit)
),
stringsAsFactors = FALSE
)
)
}
if (is.null(title) || !length(title) || is.na(title[1L]) ||
!nzchar(as.character(title)[1L])) {
title <- paste0("Meta-analysis of ", config$label)
} else {
title <- as.character(title)[1L]
}
report_dir <- tempfile("r4vn-meta-")
dir.create(report_dir, recursive = TRUE, showWarnings = FALSE)
binary_table <- .r4vn_meta_binary_table(analysis)
overview <- data.frame(
Studies = primary$k,
Effect = config$label,
Model = if (random) "Random effects" else "Fixed effect",
`Tau-squared method` = if (random) method else "Not applicable",
Inference = test_method,
stringsAsFactors = FALSE, check.names = FALSE
)
out <- list(
data = main_table,
table = main_table,
tables = list(
Overview = overview, Main = main_table,
Statistical_tests = significance_tests,
Heterogeneity = heterogeneity
),
title = title,
config = config,
effect_type = type,
analysis_data = analysis,
excluded = excluded,
model = primary,
primary_model = primary,
random_model = random_model,
fixed_model = fixed_model,
random = random,
fixed = fixed,
method = method,
hk = test_method %in% c("knha", "adhoc"),
small = small,
test_method = test_method,
small_sample = small_sample,
prediction = do_prediction,
overall = overall,
significance_tests = significance_tests,
heterogeneity = heterogeneity,
subgroup = subgroup,
subgroups = subgroup_all$analyses,
subgroup_tables = subgroup_all$wide,
subgroup_long = subgroup_all$table,
regression = regression,
bias = bias_result,
leaveout = leave_result,
influence = influence_result,
cumulative = cumulative_result,
by = by_name,
by_vars = by_names,
reg = reg_names,
cumulative_var = cumulative_name,
full = full,
profile = profile,
report = report_text,
interpretation = NULL,
overview = overview,
ci = ci,
digits = digit,
p_digits = p_digit,
plots = list(),
plot_titles = list(),
plot_display = plot_display,
plot_args = plot_args,
report_dir = report_dir,
call = match.call()
)
if (!is.null(binary_table)) out$tables$Binary_2x2 <- binary_table
if (length(subgroup_all$wide)) {
out$tables$Subgroup <- subgroup_all$wide[[1L]]
if (length(subgroup_all$wide) > 1L) {
additional_subgroups <- subgroup_all$wide[-1L]
for (nm in names(additional_subgroups)) {
out$tables[[paste0("Subgroup_", nm)]] <- additional_subgroups[[nm]]
}
}
}
if (!is.null(subgroup_all$test)) out$tables$Subgroup_test <- subgroup_all$test
if (!is.null(regression$table)) {
out$tables$Meta_regression <- regression$table
}
if (!is.null(regression$univariable)) {
out$tables$Moderator_univariable <- regression$univariable
}
if (!is.null(regression$statistics)) {
out$tables$Meta_regression_statistics <- regression$statistics
}
if (!is.null(bias_result$table)) {
out$tables$Publication_bias <- bias_result$table
}
if (!is.null(bias_result$adjusted)) {
out$tables$Publication_bias_adjusted <- bias_result$adjusted
}
if (!is.null(small_sample$comparison)) {
out$tables$Small_sample_inference <- small_sample$comparison
}
if (!is.null(leave_result$table)) {
out$tables$Leave_one_out <- leave_result$table
}
if (!is.null(influence_result$table)) {
out$tables$Influence <- influence_result$table
}
if (!is.null(cumulative_result$table)) {
out$tables$Cumulative <- cumulative_result$table
}
out$estimates <- list(
overall = overall,
subgroup = subgroup_all$table,
meta_regression = if (is.null(regression)) NULL else regression$table,
publication_bias_adjusted = if (is.null(bias_result)) NULL else bias_result$adjusted,
cumulative = if (is.null(cumulative_result)) NULL else cumulative_result$table
)
out$tests <- list(
overall = significance_tests[1L, , drop = FALSE],
heterogeneity = significance_tests[2L, , drop = FALSE],
significance = significance_tests,
subgroup = subgroup_all$test,
moderators = if (is.null(regression)) NULL else regression$statistics,
small_study_effects = if (is.null(bias_result)) NULL else bias_result$table
)
out$diagnostics <- list(
small_sample = small_sample,
publication_bias = bias_result,
leave_one_out = leave_result,
influence = influence_result
)
out$models <- list(
primary = primary,
random = random_model,
fixed = fixed_model,
subgroups = subgroup_all$analyses,
meta_regression = if (is.null(regression)) NULL else regression$model,
trimfill = if (is.null(bias_result)) NULL else bias_result$trimfill,
selection = if (is.null(bias_result)) NULL else bias_result$selection
)
out$metadata <- list(
profile = profile, effect_type = type, effect_label = config$label,
studies = primary$k, excluded = nrow(excluded), by = by_names,
moderators = reg_names, method = method, inference = test_method,
ci = ci, continuity_correction = cc, zero = zero,
variable_labels = variable_labels,
plot_display = plot_display,
input = if (all(abcd_present)) "a,b,c,d" else if (use_generic) "generic" else "raw"
)
class(out) <- c("r4vn_meta", "r4vn_tab")
out$results_text <- .r4vn_meta_results_text(out)
if (isTRUE(report_text)) {
out$interpretation <- .r4vn_meta_interpretation(out)
out$tables$Interpretation <- out$interpretation
}
if (isTRUE(plot)) {
plot_opts <- function(type) {
nested <- any(names(plot_args) %in% c(
"all", "forest", "funnel", "trimfill", "leaveout", "influence",
"baujat", "radial", "labbe", "cumulative", "bubble"
))
opts <- if (nested) {
utils::modifyList(plot_args$all %||% list(), plot_args[[type]] %||% list())
} else plot_args
opts[c("x", "type", "file")] <- NULL
opts
}
save_plot <- function(type, path, extra = list(), force = list(),
object = out) {
opts <- utils::modifyList(extra, plot_opts(type))
opts <- utils::modifyList(opts, force)
args <- c(list(x = object, type = type, file = path), opts)
ok <- tryCatch(
{
do.call(.r4vn_meta_plot_file, args)
TRUE
},
error = function(e) {
if (isTRUE(strict)) stop(e)
warning("Meta-analysis plot `", type, "` was unavailable: ",
conditionMessage(e), call. = FALSE)
FALSE
}
)
if (!isTRUE(ok) && file.exists(path)) unlink(path)
invisible(ok)
}
subgroup_plot_object <- function(subgroup_analysis, group_name) {
idx <- subgroup_analysis$indices[[group_name]]
if (is.null(idx) || sum(idx) < 2L) return(NULL)
subgroup_object <- out
subgroup_object$analysis_data <- analysis[idx, , drop = FALSE]
subgroup_object$primary_model <- subgroup_analysis$fits[[group_name]]
subgroup_object$model <- subgroup_object$primary_model
subgroup_object$random_model <- if (random) {
subgroup_object$primary_model
} else NULL
subgroup_object$fixed_model <- if (fixed && !random) {
subgroup_object$primary_model
} else NULL
subgroup_object$by <- NULL
subgroup_object$title <- paste0(
"Subgroup forest plot: ",
subgroup_analysis$label %||% subgroup_analysis$variable,
" = ", group_name
)
subgroup_object
}
out$plots$forest <- file.path(report_dir, "forest.png")
out$plot_titles$forest <- "Forest plot"
save_plot("forest", out$plots$forest)
if (length(out$subgroups)) {
subgroup_number <- 0L
for (moderator_name in names(out$subgroups)) {
subgroup_analysis <- out$subgroups[[moderator_name]]
if (!length(subgroup_analysis$fits)) next
for (group_name in names(subgroup_analysis$fits)) {
subgroup_number <- subgroup_number + 1L
subgroup_object <- subgroup_plot_object(
subgroup_analysis, group_name
)
if (is.null(subgroup_object)) next
subgroup_title <- subgroup_object$title
plot_name <- paste0("forest_subgroup_", subgroup_number)
out$plots[[plot_name]] <- file.path(
report_dir, paste0(plot_name, ".png")
)
out$plot_titles[[plot_name]] <- subgroup_title
save_plot(
"forest", out$plots[[plot_name]],
force = list(title = subgroup_title), object = subgroup_object
)
}
}
}
if (do_bias && primary$k >= 3L) {
out$plots$funnel <- file.path(report_dir, "funnel.png")
out$plot_titles$funnel <- "Contour-enhanced funnel plot"
save_plot("funnel", out$plots$funnel, list(contour = TRUE))
if (!is.null(bias_result$trimfill)) {
out$plots$trimfill <- file.path(report_dir, "trimfill.png")
out$plot_titles$trimfill <- "Trim-and-fill funnel plot"
save_plot("trimfill", out$plots$trimfill, list(contour = TRUE))
}
}
if (do_leaveout && !is.null(leave_result$table)) {
out$plots$leaveout <- file.path(report_dir, "leaveout.png")
out$plot_titles$leaveout <- "Leave-one-out analysis"
save_plot("leaveout", out$plots$leaveout)
}
if (do_influence && !is.null(influence_result$table)) {
out$plots$influence <- file.path(report_dir, "influence.png")
out$plot_titles$influence <- "Influence diagnostics"
save_plot("influence", out$plots$influence)
out$plots$baujat <- file.path(report_dir, "baujat.png")
out$plot_titles$baujat <- "Baujat plot"
save_plot("baujat", out$plots$baujat)
}
if (!is.null(cumulative_result$table)) {
out$plots$cumulative <- file.path(report_dir, "cumulative.png")
out$plot_titles$cumulative <- paste0(
"Cumulative meta-analysis: ",
cumulative_result$label %||% cumulative_name
)
save_plot(
"cumulative", out$plots$cumulative,
list(title = out$plot_titles$cumulative)
)
}
if (!is.null(regression$model) && length(reg_names) == 1L) {
moderator_label <- .r4vn_meta_variable_label(
analysis[[reg_names[1L]]], reg_names[1L]
)
out$plots$bubble <- file.path(report_dir, "bubble.png")
out$plot_titles$bubble <- paste0("Meta-regression: ", moderator_label)
save_plot("bubble", out$plots$bubble,
list(moderator = reg_names[1L],
title = out$plot_titles$bubble))
}
if (identical(profile, "full")) {
out$plots$radial <- file.path(report_dir, "radial.png")
out$plot_titles$radial <- "Radial plot"
save_plot("radial", out$plots$radial)
if (!is.null(binary_table)) {
out$plots$labbe <- file.path(report_dir, "labbe.png")
out$plot_titles$labbe <- "L'Abbe plot"
save_plot("labbe", out$plots$labbe)
}
}
available <- vapply(out$plots, file.exists, logical(1))
out$plots <- out$plots[available]
out$plot_titles <- out$plot_titles[names(out$plots)]
}
out$html <- .r4vn_meta_report_html(out)
out$table_html <- .r4vn_meta_html_table(out$table, out$title)
out$file <- file.path(report_dir, "index.html")
writeLines(enc2utf8(out$html), out$file, useBytes = TRUE)
if (isTRUE(plot) && interactive() && length(plot_display)) {
display_types <- if ("all" %in% plot_display) {
intersect(names(out$plots), plot_types)
} else plot_display
missing_display <- setdiff(display_types, names(out$plots))
if (length(missing_display)) {
warning(
"Plot-pane figure(s) unavailable for this analysis: ",
paste(missing_display, collapse = ", "), ".",
call. = FALSE
)
}
display_types <- intersect(display_types, names(out$plots))
for (display_type in display_types) {
display_args <- c(
list(x = out, type = display_type), plot_opts(display_type)
)
display_args[c("file", "width", "height", "res")] <- NULL
tryCatch(
do.call(.r4vn_meta_draw, display_args),
error = function(e) {
if (isTRUE(strict)) stop(e)
warning(
"Meta-analysis Plot-pane figure `", display_type,
"` was unavailable: ", conditionMessage(e), call. = FALSE
)
}
)
if (identical(display_type, "forest") && length(out$subgroups)) {
subgroup_drawn <- 0L
for (moderator_name in names(out$subgroups)) {
subgroup_analysis <- out$subgroups[[moderator_name]]
if (!length(subgroup_analysis$fits)) next
for (group_name in names(subgroup_analysis$fits)) {
subgroup_object <- subgroup_plot_object(
subgroup_analysis, group_name
)
if (is.null(subgroup_object)) next
subgroup_args <- c(
list(x = subgroup_object, type = "forest"),
plot_opts("forest")
)
subgroup_args[c("file", "width", "height", "res")] <- NULL
# A subgroup title must remain explicit even if a common forest
# title was supplied through plot_args.
subgroup_args$title <- subgroup_object$title
drawn <- tryCatch(
{
do.call(.r4vn_meta_draw, subgroup_args)
TRUE
},
error = function(e) {
if (isTRUE(strict)) stop(e)
warning(
"Subgroup forest plot `", subgroup_object$title,
"` was unavailable in the Plot pane: ",
conditionMessage(e), call. = FALSE
)
FALSE
}
)
if (isTRUE(drawn)) subgroup_drawn <- subgroup_drawn + 1L
}
}
if (subgroup_drawn > 0L) {
message(
"R4VN added the overall forest plot and ", subgroup_drawn,
" subgroup forest plot(s) to Plot history. Use Previous/Next ",
"to review them."
)
}
}
}
}
if (!is.null(export) || !is.null(file)) {
export_format <- export
if (is.null(export_format)) {
export_format <- tools::file_ext(as.character(file)[1L])
if (!nzchar(export_format)) export_format <- "html"
}
out$export <- .r4vn_export_meta(
out, export = export_format, file = file,
open = open, title = title
)
}
if (isTRUE(show)) .r4vn_meta_open(out$file)
invisible(out)
}
#' Print or reopen an R4VN meta-analysis
#'
#' @param x An object created by `tabmeta()`.
#' @param ... Additional arguments ignored.
#' @return The input object invisibly.
#' @export
print.r4vn_meta <- function(x, ...) {
if (!inherits(x, "r4vn_meta")) {
stop("`x` must be created by `tabmeta()`.", call. = FALSE)
}
if (!is.null(x$file) && file.exists(x$file)) {
.r4vn_meta_open(x$file)
} else {
cat(x$title, "\n")
print(x$table, row.names = FALSE)
cat("\n")
print(x$significance_tests, row.names = FALSE)
cat("\n")
print(x$heterogeneity, row.names = FALSE)
}
invisible(x)
}
#' Convert an R4VN meta-analysis to a data frame
#'
#' @param x An object created by `tabmeta()`.
#' @param row.names Ignored.
#' @param optional Ignored.
#' @param ... Additional arguments ignored.
#' @return The main publication table.
#' @export
as.data.frame.r4vn_meta <- function(x, row.names = NULL,
optional = FALSE, ...) {
x$table
}
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.