R/tabforest.R

In R4VN: Health Data Analysis and Publication-Ready Reporting

# =============================================================================
# R4VN tabforest() - regression forest table/plot engine
# =============================================================================

.r4vn_tf_escape_name <- function(x) {
  x <- as.character(x)
  ifelse(make.names(x) == x, x, paste0("`", gsub("`", "", x, fixed = TRUE), "`"))
}

.r4vn_tf_data <- function(data = NULL) {
  if (!is.null(data)) {
    if (!is.data.frame(data)) stop("`data` must be a data frame.", call. = FALSE)
    return(data)
  }

  # Use the same active-data resolver as the current R4VN analysis commands.
  if (exists(".r4vn_resolve_analysis_data", mode = "function", inherits = TRUE)) {
    z <- tryCatch(
      get(".r4vn_resolve_analysis_data", mode = "function", inherits = TRUE)(NULL),
      error = function(e) NULL
    )
    if (is.data.frame(z)) return(z)
  }
  if (exists(".r4vn_get_active", mode = "function", inherits = TRUE)) {
    f <- get(".r4vn_get_active", mode = "function", inherits = TRUE)
    z <- tryCatch(f(required = FALSE), error = function(e) tryCatch(f(), error = function(e2) NULL))
    if (is.data.frame(z)) return(z)
  }

  # Backward-compatible fallbacks for older development builds.
  if (exists("active_data", mode = "function", inherits = TRUE)) {
    f <- get("active_data", mode = "function", inherits = TRUE)
    z <- tryCatch(f(), error = function(e) NULL)
    if (is.data.frame(z)) return(z)
  }
  nm <- getOption(".r4vn_active_data_name")
  if (is.character(nm) && length(nm) == 1L && nzchar(nm)) {
    z <- tryCatch(get(nm, envir = .GlobalEnv, inherits = TRUE), error = function(e) NULL)
    if (is.data.frame(z)) return(z)
  }
  stop("No data frame was supplied and no active R4VN data frame is available. Use `data = ...` or `usedf()` first.", call. = FALSE)
}

.r4vn_tf_name <- function(expr, data, env, arg, optional = FALSE) {
  if (identical(expr, quote(NULL))) {
    if (optional) return(NULL)
    stop(sprintf("`%s` is required.", arg), call. = FALSE)
  }
  if (is.symbol(expr)) {
    nm <- as.character(expr)
    if (nm %in% names(data)) return(nm)
  }
  value <- tryCatch(eval(expr, envir = env), error = function(e) NULL)
  if (is.character(value) && length(value) == 1L && value %in% names(data)) return(value)
  raw <- paste(deparse(expr, width.cutoff = 500L), collapse = "")
  if (raw %in% names(data)) return(raw)
  stop(sprintf("`%s` must identify one variable in `data`.", arg), call. = FALSE)
}

.r4vn_tf_levels <- function(x) {
  observed <- x[!is.na(x)]
  if (is.factor(x)) {
    lev <- levels(x)
    return(lev[lev %in% as.character(observed)])
  }
  if (is.logical(x)) return(as.character(c(FALSE, TRUE)[c(FALSE, TRUE) %in% observed]))
  z <- unique(as.character(observed))
  if (is.numeric(x)) z <- as.character(sort(unique(x[!is.na(x)])))
  z
}

.r4vn_tf_spec <- function(x, data, arg = "predictors", allow_null = FALSE) {
  if (is.null(x) || identical(x, FALSE)) {
    if (allow_null) return(NULL)
    stop(sprintf("`%s` is required.", arg), call. = FALSE)
  }
  if (inherits(x, "r4vn_vars")) {
    resolver <- get0(".r4vn_resolve_vars", mode = "function", inherits = TRUE)
    out <- if (is.function(resolver)) {
      as.data.frame(resolver(x, data = data, default_type = "auto", strict = TRUE), stringsAsFactors = FALSE)
    } else {
      as.data.frame(x, stringsAsFactors = FALSE)
    }
  } else if (is.character(x)) {
    out <- lapply(x, function(v) {
      if (!v %in% names(data)) stop(sprintf("Variable `%s` in `%s` was not found in `data`.", v, arg), call. = FALSE)
      is_cont <- is.numeric(data[[v]])
      data.frame(variable = v,
                 type = if (is_cont) "mean" else "categorical",
                 specification = v,
                 reference_index = if (is_cont) NA_integer_ else 1L,
                 stringsAsFactors = FALSE)
    })
    out <- do.call(rbind, out)
  } else {
    stop(sprintf("`%s` must be created by `vars()` or supplied as a character vector of variable names.", arg), call. = FALSE)
  }
  need <- c("variable", "type", "specification", "reference_index")
  if (!all(need %in% names(out))) stop(sprintf("`%s` is not a valid R4VN variable specification.", arg), call. = FALSE)
  if (any(!out$variable %in% names(data))) {
    stop(sprintf("Variable(s) not found in `data`: %s.", paste(out$variable[!out$variable %in% names(data)], collapse = ", ")), call. = FALSE)
  }
  out <- out[!duplicated(out$variable), need, drop = FALSE]
  rownames(out) <- NULL
  out
}

.r4vn_tf_merge_spec <- function(...) {
  z <- Filter(function(x) !is.null(x) && nrow(x), list(...))
  if (!length(z)) return(NULL)
  out <- do.call(rbind, z)
  out <- out[!duplicated(out$variable, fromLast = TRUE), , drop = FALSE]
  rownames(out) <- NULL
  out
}

.r4vn_tf_ref <- function(x, spec_row) {
  lev <- .r4vn_tf_levels(x)
  idx <- suppressWarnings(as.integer(spec_row$reference_index[1L]))
  if (!length(lev)) return(NULL)
  if (is.na(idx)) idx <- 1L
  if (idx < 1L || idx > length(lev)) {
    stop(sprintf("Reference level b%s is invalid for `%s`, which has %s observed level(s).",
                 idx, spec_row$variable[1L], length(lev)), call. = FALSE)
  }
  lev[idx]
}

.r4vn_tf_var_label <- function(data, variable, labels = NULL) {
  if (!is.null(labels) && !is.null(names(labels)) && variable %in% names(labels)) {
    z <- as.character(labels[[variable]])[1L]
    if (nzchar(z)) return(z)
  }
  z <- attr(data[[variable]], "label", exact = TRUE)
  if (!is.null(z) && length(z) && nzchar(as.character(z)[1L])) return(as.character(z)[1L])
  variable
}

.r4vn_tf_level_label <- function(variable, level, level_labels = NULL) {
  if (!is.null(level_labels) && !is.null(level_labels[[variable]])) {
    z <- level_labels[[variable]]
    if (!is.null(names(z)) && level %in% names(z)) return(as.character(z[[level]])[1L])
  }
  as.character(level)
}

.r4vn_tf_per <- function(variable, per = NULL) {
  if (is.null(per)) return(1)
  if (length(per) == 1L && is.null(names(per))) return(as.numeric(per)[1L])
  if (!is.null(names(per)) && variable %in% names(per)) return(as.numeric(per[[variable]])[1L])
  1
}

.r4vn_tf_per_label <- function(variable, value, per_labels = NULL, lang = "en") {
  if (!is.null(per_labels) && !is.null(names(per_labels)) && variable %in% names(per_labels)) {
    return(as.character(per_labels[[variable]])[1L])
  }
  if (!is.finite(value) || value == 1) return(NULL)
  if (identical(lang, "vi")) paste0("m\u1ed7i ", format(value, trim = TRUE, scientific = FALSE), " \u0111\u01a1n v\u1ecb")
  else paste0("per ", format(value, trim = TRUE, scientific = FALSE), " units")
}

.r4vn_tf_prepare <- function(data, outcome_name, time_name, status_name,
                             model_spec, effect, event = NULL, failure = NULL) {
  vars <- if (is.null(model_spec)) character() else model_spec$variable
  needed <- unique(c(outcome_name, time_name, status_name, vars))
  needed <- needed[!is.na(needed) & nzchar(needed)]
  d <- data[, needed, drop = FALSE]
  keep <- stats::complete.cases(d)
  d <- d[keep, , drop = FALSE]
  if (!nrow(d)) stop("No complete observations are available for this model.", call. = FALSE)

  if (identical(effect, "HR")) {
    d$.time <- suppressWarnings(as.numeric(d[[time_name]]))
    status <- d[[status_name]]
    lev <- .r4vn_tf_levels(status)
    if (is.null(failure)) {
      if (is.numeric(status) && any(status == 1, na.rm = TRUE)) failure <- 1
      else failure <- tail(lev, 1L)
    }
    d$.status <- as.integer(as.character(status) == as.character(failure))
    if (any(!is.finite(d$.time)) || any(d$.time < 0)) stop("Survival time must be non-negative and finite.", call. = FALSE)
    if (length(unique(d$.status)) < 2L) stop("The survival status contains fewer than two event states in the analysis sample.", call. = FALSE)
  } else if (effect %in% c("OR", "RR", "PR")) {
    y <- d[[outcome_name]]
    lev <- .r4vn_tf_levels(y)
    if (length(lev) != 2L) stop(sprintf("%s requires a binary outcome.", effect), call. = FALSE)
    if (is.null(event)) event <- tail(lev, 1L)
    if (!as.character(event) %in% lev) stop("`event` is not an observed level of the outcome.", call. = FALSE)
    d$.outcome <- as.integer(as.character(y) == as.character(event))
  } else if (identical(effect, "IRR")) {
    y <- suppressWarnings(as.numeric(d[[outcome_name]]))
    if (any(!is.finite(y)) || any(y < 0)) stop("IRR requires a non-negative count outcome.", call. = FALSE)
    if (any(abs(y - round(y)) > sqrt(.Machine$double.eps))) stop("IRR requires integer count values.", call. = FALSE)
    d$.outcome <- y
  } else {
    y <- suppressWarnings(as.numeric(d[[outcome_name]]))
    if (any(!is.finite(y))) stop("Beta regression requires a numeric continuous outcome.", call. = FALSE)
    d$.outcome <- y
  }

  if (!is.null(model_spec) && nrow(model_spec)) {
    for (i in seq_len(nrow(model_spec))) {
      v <- model_spec$variable[i]
      if (identical(model_spec$type[i], "categorical")) {
        lev <- .r4vn_tf_levels(data[[v]])
        ref <- .r4vn_tf_ref(data[[v]], model_spec[i, , drop = FALSE])
        d[[v]] <- factor(d[[v]], levels = lev)
        d[[v]] <- stats::relevel(d[[v]], ref = ref)
        if (nlevels(droplevels(d[[v]])) < 2L) stop(sprintf("Categorical variable `%s` has fewer than two observed levels in this model.", v), call. = FALSE)
      } else {
        d[[v]] <- suppressWarnings(as.numeric(d[[v]]))
        if (!is.finite(stats::sd(d[[v]])) || stats::sd(d[[v]]) == 0) stop(sprintf("Continuous variable `%s` has no variation in this model.", v), call. = FALSE)
      }
    }
  }
  attr(d, "r4vn_keep") <- which(keep)
  attr(d, "r4vn_event") <- event
  attr(d, "r4vn_failure") <- failure
  d
}

.r4vn_tf_robust_poisson <- function(fit) {
  X <- stats::model.matrix(fit)
  mu <- stats::fitted(fit)
  y <- fit$y
  if (is.null(y)) y <- tryCatch(stats::model.response(stats::model.frame(fit)), error = function(e) NULL)
  if (is.null(y)) stop("The Poisson model does not retain a usable response.", call. = FALSE)
  r <- y - mu
  bread <- tryCatch(solve(crossprod(X, X * as.vector(mu))), error = function(e) NULL)
  if (is.null(bread)) bread <- tryCatch(qr.solve(crossprod(X, X * as.vector(mu))), error = function(e) NULL)
  if (is.null(bread)) stop("The robust covariance matrix could not be estimated.", call. = FALSE)
  meat <- crossprod(X, X * as.vector(r^2))
  V <- bread %*% meat %*% bread
  dimnames(V) <- list(colnames(X), colnames(X))
  V
}

.r4vn_tf_formula <- function(response, variables) {
  terms <- vapply(variables, .r4vn_tf_escape_name, character(1))
  stats::as.formula(paste(response, "~", if (length(terms)) paste(terms, collapse = " + ") else "1"))
}

.r4vn_tf_fit <- function(d, model_spec, effect, ci = .95) {
  variables <- if (is.null(model_spec)) character() else model_spec$variable
  if (identical(effect, "HR")) {
    if (!requireNamespace("survival", quietly = TRUE)) {
      stop("Cox forest plots require the `survival` package. Install it with install.packages('survival').", call. = FALSE)
    }
    f <- .r4vn_tf_formula("survival::Surv(.time, .status)", variables)
    fit <- survival::coxph(f, data = d, x = TRUE, y = TRUE, model = TRUE, ties = "efron")
    V <- stats::vcov(fit)
    family <- "cox"
  } else if (identical(effect, "Beta")) {
    f <- .r4vn_tf_formula(".outcome", variables)
    fit <- stats::lm(f, data = d, x = TRUE, y = TRUE)
    V <- stats::vcov(fit)
    family <- "lm"
  } else if (identical(effect, "OR")) {
    f <- .r4vn_tf_formula(".outcome", variables)
    fit <- stats::glm(f, family = stats::binomial("logit"), data = d, x = TRUE, y = TRUE)
    V <- stats::vcov(fit)
    family <- "glm"
  } else {
    f <- .r4vn_tf_formula(".outcome", variables)
    fit <- stats::glm(f, family = stats::poisson("log"), data = d, x = TRUE, y = TRUE)
    V <- if (effect %in% c("RR", "PR")) .r4vn_tf_robust_poisson(fit) else stats::vcov(fit)
    family <- if (effect %in% c("RR", "PR")) "poisson_robust" else "poisson"
  }
  list(fit = fit, vcov = V, family = family, data = d, spec = model_spec, ci = ci)
}

.r4vn_tf_term_index <- function(fit, variable) {
  tt <- attr(stats::terms(fit), "term.labels")
  clean <- gsub("`", "", tt, fixed = TRUE)
  match(variable, clean)
}

.r4vn_tf_coef_names <- function(fit, variable) {
  idx <- .r4vn_tf_term_index(fit, variable)
  if (is.na(idx)) return(character())
  mm <- stats::model.matrix(fit)
  ass <- attr(mm, "assign")
  colnames(mm)[ass == idx]
}

.r4vn_tf_global_p <- function(model, variable) {
  b <- stats::coef(model$fit)
  cn <- intersect(.r4vn_tf_coef_names(model$fit, variable), names(b))
  cn <- cn[is.finite(b[cn])]
  if (!length(cn)) return(NA_real_)
  V <- model$vcov[cn, cn, drop = FALSE]
  good <- is.finite(b[cn]) & is.finite(diag(V)) & diag(V) > 0
  cn <- cn[good]
  if (!length(cn)) return(NA_real_)
  b <- b[cn]
  V <- V[cn, cn, drop = FALSE]
  inv <- tryCatch(solve(V), error = function(e) tryCatch(qr.solve(V), error = function(e2) NULL))
  if (is.null(inv)) return(NA_real_)
  q <- as.numeric(t(b) %*% inv %*% b)
  df <- qr(V)$rank
  if (!is.finite(q) || df < 1L) return(NA_real_)
  stats::pchisq(q, df = df, lower.tail = FALSE)
}

.r4vn_tf_one_effect <- function(model, variable, spec_row, effect, per = 1,
                                ci = .95, model_name = "Crude") {
  fit <- model$fit
  V <- model$vcov
  b <- stats::coef(fit)
  cn <- intersect(.r4vn_tf_coef_names(fit, variable), names(b))
  alpha <- 1 - ci
  crit <- if (inherits(fit, "lm")) stats::qt(1 - alpha / 2, df = stats::df.residual(fit)) else stats::qnorm(1 - alpha / 2)
  n <- stats::nobs(fit)
  events <- if (identical(effect, "HR") && ".status" %in% names(model$data)) {
    sum(model$data$.status == 1L)
  } else if (effect %in% c("OR", "RR", "PR") && ".outcome" %in% names(model$data)) {
    sum(model$data$.outcome == 1L)
  } else NA_real_
  gp <- .r4vn_tf_global_p(model, variable)

  if (identical(spec_row$type[1L], "categorical")) {
    lev <- levels(model$data[[variable]])
    ref <- lev[1L]
    nonref <- lev[-1L]
    out <- data.frame(variable = variable, level = lev, reference = lev == ref,
                      model = model_name, estimate = NA_real_, lower = NA_real_, upper = NA_real_,
                      p = NA_real_, global_p = gp, n = n, events = events,
                      stringsAsFactors = FALSE)
    if (effect %in% c("OR", "RR", "PR", "IRR", "HR")) out$estimate[out$reference] <- 1
    usable <- cn[is.finite(b[cn])]
    k <- min(length(usable), length(nonref))
    if (k) {
      for (j in seq_len(k)) {
        nm <- usable[j]
        beta <- unname(b[nm])
        se <- sqrt(unname(V[nm, nm]))
        if (!is.finite(se) || se <= 0) next
        z <- beta / se
        p <- if (inherits(fit, "lm")) 2 * stats::pt(abs(z), df = stats::df.residual(fit), lower.tail = FALSE) else 2 * stats::pnorm(abs(z), lower.tail = FALSE)
        lo <- beta - crit * se
        hi <- beta + crit * se
        row <- which(out$level == nonref[j])[1L]
        if (identical(effect, "Beta")) {
          out$estimate[row] <- beta
          out$lower[row] <- lo
          out$upper[row] <- hi
        } else {
          out$estimate[row] <- exp(beta)
          out$lower[row] <- exp(lo)
          out$upper[row] <- exp(hi)
        }
        out$p[row] <- p
      }
    }
    return(out)
  }

  if (!length(cn)) return(NULL)
  nm <- cn[1L]
  beta <- unname(b[nm])
  se <- sqrt(unname(V[nm, nm]))
  if (!is.finite(beta) || !is.finite(se) || se <= 0) return(NULL)
  mult <- if (is.finite(per) && per > 0) per else 1
  z <- beta / se
  p <- if (inherits(fit, "lm")) 2 * stats::pt(abs(z), df = stats::df.residual(fit), lower.tail = FALSE) else 2 * stats::pnorm(abs(z), lower.tail = FALSE)
  lo <- beta - crit * se
  hi <- beta + crit * se
  if (identical(effect, "Beta")) {
    est <- beta * mult; lower <- lo * mult; upper <- hi * mult
  } else {
    est <- exp(beta * mult); lower <- exp(lo * mult); upper <- exp(hi * mult)
  }
  data.frame(variable = variable, level = "", reference = FALSE,
             model = model_name, estimate = est, lower = lower, upper = upper,
             p = p, global_p = gp, n = n, events = events,
             stringsAsFactors = FALSE)
}

.r4vn_tf_effect_type <- function(data, outcome_name, time_name,
                                 or = FALSE, rr = FALSE, pr = FALSE, irr = FALSE,
                                 estimate = c("auto", "beta", "or", "rr", "pr", "irr", "hr")) {
  estimate <- match.arg(tolower(estimate[1L]), c("auto", "beta", "or", "rr", "pr", "irr", "hr"))
  flags <- c(or = isTRUE(or), rr = isTRUE(rr), pr = isTRUE(pr), irr = isTRUE(irr))
  if (sum(flags) > 1L) stop("Choose only one of `or`, `rr`, `pr`, or `irr`.", call. = FALSE)
  if (!is.null(time_name)) {
    if (estimate != "auto" && estimate != "hr") stop("When `time` is supplied, `estimate` must be 'auto' or 'hr'.", call. = FALSE)
    if (any(flags)) stop("Do not use `or`, `rr`, `pr`, or `irr` with a survival outcome.", call. = FALSE)
    return("HR")
  }
  if (estimate != "auto") return(switch(estimate, beta = "Beta", or = "OR", rr = "RR", pr = "PR", irr = "IRR", hr = "HR"))
  if (any(flags)) return(unname(c(OR = "OR", RR = "RR", PR = "PR", IRR = "IRR")[toupper(names(flags)[which(flags)])][1L]))
  y <- data[[outcome_name]]
  if (length(.r4vn_tf_levels(y)) == 2L) return("OR")
  if (is.numeric(y)) return("Beta")
  stop("The outcome is neither binary nor numeric. Specify a supported `estimate` or recode the outcome.", call. = FALSE)
}

.r4vn_tf_text <- function(lang = c("en", "vi"), text = NULL) {
  lang <- match.arg(lang)
  out <- if (lang == "vi") list(
    characteristic = "Y\u1ebfu t\u1ed1",
    reference = "Tham chi\u1ebfu",
    crude = "\u0110\u01a1n bi\u1ebfn",
    adjusted = "Hi\u1ec7u ch\u1ec9nh",
    multi = "\u0110a bi\u1ebfn",
    p = "p",
    n = "n",
    events = "Bi\u1ebfn c\u1ed1",
    arrow_note = "M\u0169i t\u00ean cho bi\u1ebft KTC v\u01b0\u1ee3t ra ngo\u00e0i gi\u1edbi h\u1ea1n tr\u1ee5c.",
    beta = "Beta (KTC 95%)",
    OR = "OR (KTC 95%)",
    RR = "RR (KTC 95%)",
    PR = "PR (KTC 95%)",
    IRR = "IRR (KTC 95%)",
    HR = "HR (KTC 95%)"
  ) else list(
    characteristic = "Characteristic",
    reference = "Reference",
    crude = "Crude",
    adjusted = "Adjusted",
    multi = "Multivariable",
    p = "p",
    n = "n",
    events = "Events",
    arrow_note = "Arrows indicate confidence intervals extending beyond the plotting range.",
    beta = "Beta (95% CI)",
    OR = "OR (95% CI)",
    RR = "RR (95% CI)",
    PR = "PR (95% CI)",
    IRR = "IRR (95% CI)",
    HR = "HR (95% CI)"
  )
  if (!is.null(text)) {
    if (!is.list(text)) stop("`text` must be a named list.", call. = FALSE)
    for (nm in names(text)) out[[nm]] <- as.character(text[[nm]])[1L]
  }
  out
}

.r4vn_tf_rows <- function(data, focal_spec, estimates, labels = NULL,
                          level_labels = NULL, per = NULL, per_labels = NULL,
                          lang = "en", reference = TRUE) {
  rows <- list()
  for (i in seq_len(nrow(focal_spec))) {
    v <- focal_spec$variable[i]
    vl <- .r4vn_tf_var_label(data, v, labels)
    if (identical(focal_spec$type[i], "categorical")) {
      lev <- .r4vn_tf_levels(data[[v]])
      rows[[length(rows) + 1L]] <- data.frame(variable = v, level = NA_character_,
        row_type = "header", label = vl, stringsAsFactors = FALSE)
      ref0 <- .r4vn_tf_ref(data[[v]], focal_spec[i, , drop = FALSE])
      for (z in lev) {
        if (!isTRUE(reference) && identical(as.character(z), as.character(ref0))) next
        rows[[length(rows) + 1L]] <- data.frame(variable = v, level = as.character(z),
          row_type = "level", label = .r4vn_tf_level_label(v, as.character(z), level_labels), stringsAsFactors = FALSE)
      }
    } else {
      pv <- .r4vn_tf_per(v, per)
      pl <- .r4vn_tf_per_label(v, pv, per_labels, lang)
      if (!is.null(pl) && nzchar(pl)) vl <- paste0(vl, " (", pl, ")")
      rows[[length(rows) + 1L]] <- data.frame(variable = v, level = "",
        row_type = "continuous", label = vl, stringsAsFactors = FALSE)
    }
  }
  out <- do.call(rbind, rows)
  rownames(out) <- NULL
  out
}

.r4vn_tf_build_raw <- function(data, outcome_name, time_name, focal_spec,
                               adjusted, multi, effect, event, failure, ci,
                               crude, sample, per) {
  adjusted_spec <- NULL
  adjusted_all <- isTRUE(adjusted)
  if (!isFALSE(adjusted) && !is.null(adjusted) && !isTRUE(adjusted)) adjusted_spec <- .r4vn_tf_spec(adjusted, data, "adjusted")
  multi_spec <- NULL
  if (!isFALSE(multi) && !is.null(multi) && !isTRUE(multi)) multi_spec <- .r4vn_tf_spec(multi, data, "multi")
  if (isTRUE(multi)) multi_spec <- focal_spec

  n_model_groups <- sum(c(isTRUE(crude), !isFALSE(adjusted) && !is.null(adjusted), !isFALSE(multi) && !is.null(multi)))
  sample <- match.arg(sample, c("auto", "common", "model"))
  common <- identical(sample, "common") || (identical(sample, "auto") && n_model_groups > 1L)

  union_spec <- focal_spec
  if (adjusted_all) union_spec <- .r4vn_tf_merge_spec(union_spec, focal_spec)
  else union_spec <- .r4vn_tf_merge_spec(union_spec, adjusted_spec)
  union_spec <- .r4vn_tf_merge_spec(union_spec, multi_spec)
  analysis_data <- data
  if (common) {
    needed <- unique(c(outcome_name, time_name, if (identical(effect, "HR")) outcome_name else NULL,
                       union_spec$variable))
    needed <- needed[!is.na(needed) & nzchar(needed)]
    keep <- stats::complete.cases(data[, needed, drop = FALSE])
    analysis_data <- data[keep, , drop = FALSE]
    if (!nrow(analysis_data)) stop("`sample = 'common'` leaves no complete observations.", call. = FALSE)
  }

  estimates <- list(); models <- list(); model_keys <- character()
  status_name <- if (identical(effect, "HR")) outcome_name else NULL

  if (isTRUE(crude)) {
    model_keys <- c(model_keys, "Crude")
    models$Crude <- list()
    for (i in seq_len(nrow(focal_spec))) {
      fs <- focal_spec[i, , drop = FALSE]
      d <- .r4vn_tf_prepare(analysis_data, outcome_name, time_name, status_name, fs, effect, event, failure)
      m <- .r4vn_tf_fit(d, fs, effect, ci)
      z <- .r4vn_tf_one_effect(m, fs$variable, fs, effect, .r4vn_tf_per(fs$variable, per), ci, "Crude")
      if (!is.null(z)) estimates[[length(estimates) + 1L]] <- z
      models$Crude[[fs$variable]] <- m$fit
    }
  }

  if (!isFALSE(adjusted) && !is.null(adjusted)) {
    model_keys <- c(model_keys, "Adjusted")
    models$Adjusted <- list()
    for (i in seq_len(nrow(focal_spec))) {
      fs <- focal_spec[i, , drop = FALSE]
      adj <- if (adjusted_all) focal_spec[focal_spec$variable != fs$variable, , drop = FALSE] else adjusted_spec
      if (!is.null(adj) && nrow(adj)) adj <- adj[adj$variable != fs$variable, , drop = FALSE]
      ms <- .r4vn_tf_merge_spec(fs, adj)
      d <- .r4vn_tf_prepare(analysis_data, outcome_name, time_name, status_name, ms, effect, event, failure)
      m <- .r4vn_tf_fit(d, ms, effect, ci)
      z <- .r4vn_tf_one_effect(m, fs$variable, fs, effect, .r4vn_tf_per(fs$variable, per), ci, "Adjusted")
      if (!is.null(z)) estimates[[length(estimates) + 1L]] <- z
      models$Adjusted[[fs$variable]] <- m$fit
    }
  }

  if (!isFALSE(multi) && !is.null(multi)) {
    model_keys <- c(model_keys, "Multivariable")
    ms <- multi_spec
    d <- .r4vn_tf_prepare(analysis_data, outcome_name, time_name, status_name, ms, effect, event, failure)
    m <- .r4vn_tf_fit(d, ms, effect, ci)
    models$Multivariable <- m$fit
    for (i in seq_len(nrow(focal_spec))) {
      fs <- focal_spec[i, , drop = FALSE]
      if (!fs$variable %in% ms$variable) next
      fit_spec <- ms[match(fs$variable, ms$variable), , drop = FALSE]
      z <- .r4vn_tf_one_effect(m, fs$variable, fit_spec, effect, .r4vn_tf_per(fs$variable, per), ci, "Multivariable")
      if (!is.null(z)) estimates[[length(estimates) + 1L]] <- z
    }
  }

  if (!length(estimates)) stop("No estimable effects were produced.", call. = FALSE)
  estimates <- do.call(rbind, estimates)
  rownames(estimates) <- NULL
  list(estimates = estimates, models = models, model_keys = unique(model_keys),
       common_sample = common, analysis_n = nrow(analysis_data),
       adjusted_spec = adjusted_spec, adjusted_all = adjusted_all, multi_spec = multi_spec)
}

.r4vn_tf_fit_spec <- function(fit) {
  mf <- tryCatch(stats::model.frame(fit), error = function(e) NULL)
  tt <- attr(stats::terms(fit), "term.labels")
  tt <- gsub("`", "", tt, fixed = TRUE)
  if (!length(tt)) return(NULL)
  lapply_out <- lapply(tt, function(v) {
    x <- if (!is.null(mf) && v %in% names(mf)) mf[[v]] else NULL
    categorical <- !is.null(x) && (is.factor(x) || is.character(x) || is.logical(x))
    data.frame(variable = v, type = if (categorical) "categorical" else "mean",
               specification = v, reference_index = if (categorical) 1L else NA_integer_, stringsAsFactors = FALSE)
  })
  do.call(rbind, lapply_out)
}

.r4vn_tf_from_model <- function(fit, model_name = "Model", effect = NULL, vcov_override = NULL, ci = .95) {
  if (is.null(effect)) {
    if (inherits(fit, "coxph")) effect <- "HR"
    else if (inherits(fit, "lm") && !inherits(fit, "glm")) effect <- "Beta"
    else if (inherits(fit, "glm") && identical(fit$family$family, "binomial")) effect <- "OR"
    else if (inherits(fit, "glm") && identical(fit$family$family, "poisson")) effect <- "IRR"
    else stop("Could not infer the effect type from this model.", call. = FALSE)
  }
  spec <- .r4vn_tf_fit_spec(fit)
  if (is.null(spec) || !nrow(spec)) stop("The model has no predictor terms to display.", call. = FALSE)
  d <- tryCatch(stats::model.frame(fit), error = function(e) NULL)
  if (is.null(d)) stop("The fitted model does not retain a usable model frame.", call. = FALSE)
  V <- if (!is.null(vcov_override)) vcov_override else if (effect %in% c("RR", "PR")) .r4vn_tf_robust_poisson(fit) else stats::vcov(fit)
  model <- list(fit = fit, vcov = V, data = d, ci = ci)
  est <- list()
  for (i in seq_len(nrow(spec))) {
    v <- spec$variable[i]
    z <- .r4vn_tf_one_effect(model, v, spec[i, , drop = FALSE], effect, 1, ci, model_name)
    if (!is.null(z)) est[[length(est) + 1L]] <- z
  }
  list(estimates = do.call(rbind, est), spec = spec, effect = effect, models = setNames(list(fit), model_name), model_keys = model_name, data = d)
}

.r4vn_tf_rows_from_estimates <- function(estimates, labels = NULL, level_labels = NULL, reference = TRUE) {
  vars <- unique(estimates$variable)
  out <- list()
  for (v in vars) {
    z <- estimates[estimates$variable == v, , drop = FALSE]
    levels0 <- unique(z$level[!is.na(z$level) & nzchar(z$level)])
    has_ref <- any(z$reference %in% TRUE, na.rm = TRUE)
    vl <- if (!is.null(labels) && !is.null(names(labels)) && v %in% names(labels)) as.character(labels[[v]])[1L] else v
    if (length(levels0) > 1L || has_ref) {
      out[[length(out) + 1L]] <- data.frame(variable = v, level = NA_character_, row_type = "header", label = vl, stringsAsFactors = FALSE)
      for (lv in levels0) {
        is_ref <- any(z$level == lv & z$reference %in% TRUE, na.rm = TRUE)
        if (!isTRUE(reference) && is_ref) next
        out[[length(out) + 1L]] <- data.frame(variable = v, level = lv, row_type = "level", label = .r4vn_tf_level_label(v, lv, level_labels), stringsAsFactors = FALSE)
      }
    } else {
      out[[length(out) + 1L]] <- data.frame(variable = v, level = if (length(levels0)) levels0[1L] else "", row_type = "continuous", label = vl, stringsAsFactors = FALSE)
    }
  }
  do.call(rbind, out)
}

.r4vn_tf_object <- function(x, select = NULL, ci = .95) {
  if (inherits(x, "r4vn_surv")) {
    tabs <- x$cox
    if (is.null(tabs)) stop("This `r4vn_surv` object has no Cox regression results.", call. = FALSE)
    available <- c("crude", "adjusted", "multi")
    available <- available[vapply(available, function(k) !is.null(tabs[[k]]) && nrow(tabs[[k]]) > 0L, logical(1))]
    if (!length(available)) stop("This `r4vn_surv` object has no Cox estimates to plot.", call. = FALSE)
    if (is.null(select)) select <- available
    select <- intersect(tolower(select), available)
    if (!length(select)) stop("`select` did not match available Cox components.", call. = FALSE)
    nm_map <- c(crude = "Crude", adjusted = "Adjusted", multi = "Multivariable")
    est <- list()
    for (k in select) {
      z <- tabs[[k]]
      need <- c("variable", "level", "estimate", "lower", "upper", "p")
      if (!all(need %in% names(z))) stop("The Cox result table does not contain the fields required by `tabforest()`.", call. = FALSE)
      ref <- if ("reference" %in% names(z)) as.logical(z$reference) else is.na(z$estimate)
      est[[length(est) + 1L]] <- data.frame(variable = as.character(z$variable), level = as.character(z$level),
        reference = ref, model = unname(nm_map[k]), estimate = as.numeric(z$estimate), lower = as.numeric(z$lower), upper = as.numeric(z$upper),
        p = as.numeric(z$p), global_p = if ("global_p" %in% names(z)) as.numeric(z$global_p) else NA_real_,
        n = if ("n" %in% names(z)) as.numeric(z$n) else NA_real_, events = if ("events" %in% names(z)) as.numeric(z$events) else NA_real_, stringsAsFactors = FALSE)
    }
    return(list(estimates = do.call(rbind, est), effect = "HR", models = tabs, model_keys = unname(nm_map[select]), data = NULL, spec = NULL))
  }

  if (inherits(x, "r4vn_tabmulti")) {
    fits <- x$models
    if (is.null(fits) || !length(fits)) stop("This `r4vn_tabmulti` object has no fitted models.", call. = FALSE)
    available <- names(fits)
    if (is.null(select)) select <- if ("full" %in% available) "full" else available[1L]
    select <- intersect(select, available)
    if (!length(select)) stop(sprintf("Available `tabmulti()` models are: %s.", paste(available, collapse = ", ")), call. = FALSE)
    effect <- x$effect
    if (is.null(effect) || !nzchar(as.character(effect)[1L])) {
      fit0 <- fits[[select[1L]]]
      effect <- if (inherits(fit0, "glm") && identical(fit0$family$family, "binomial")) "OR" else if (inherits(fit0, "glm") && identical(fit0$family$family, "poisson")) "RR" else "Beta"
    }
    effect <- toupper(as.character(effect)[1L])
    if (effect %in% c("BETA", "COEFFICIENT")) effect <- "Beta"
    spec <- x$metadata
    if (is.null(spec) || !nrow(spec)) spec <- .r4vn_tf_fit_spec(fits[[select[1L]]])
    est <- list(); model_list <- list()
    for (k in select) {
      fit <- fits[[k]]
      model_list[[k]] <- fit
      V <- if (effect %in% c("RR", "PR")) .r4vn_tf_robust_poisson(fit) else stats::vcov(fit)
      model <- list(fit = fit, vcov = V, data = stats::model.frame(fit), ci = ci)
      selected_vars <- if (!is.null(x$selected[[k]])) x$selected[[k]] else spec$variable
      for (i in seq_len(nrow(spec))) {
        if (!spec$variable[i] %in% selected_vars) next
        z <- .r4vn_tf_one_effect(model, spec$variable[i], spec[i, , drop = FALSE], effect, 1, ci, k)
        if (!is.null(z)) est[[length(est) + 1L]] <- z
      }
    }
    return(list(estimates = do.call(rbind, est), effect = effect, models = model_list, model_keys = select, data = NULL, spec = spec))
  }

  if (inherits(x, c("r4vn_stat", "r4vn_result")) && !is.null(x$raw$model)) {
    fit <- x$raw$model
    V <- if (!is.null(x$raw$vcov)) x$raw$vcov else NULL
    effect <- NULL
    ttl <- if (!is.null(x$title)) tolower(as.character(x$title)[1L]) else ""
    if (inherits(fit, "glm") && identical(fit$family$family, "binomial")) effect <- "OR"
    else if (inherits(fit, "glm") && identical(fit$family$family, "poisson")) effect <- "IRR"
    else if (inherits(fit, "lm") && !inherits(fit, "glm")) effect <- "Beta"
    else if (inherits(fit, "coxph")) effect <- "HR"
    if (grepl("logistic", ttl, fixed = TRUE)) effect <- "OR"
    if (grepl("poisson", ttl, fixed = TRUE)) effect <- "IRR"
    if (grepl("linear regression", ttl, fixed = TRUE)) effect <- "Beta"
    return(.r4vn_tf_from_model(fit, "Model", effect, V, ci))
  }

  if (inherits(x, c("lm", "glm", "coxph"))) return(.r4vn_tf_from_model(x, "Model", NULL, NULL, ci))
  stop("Unsupported object. `tabforest()` currently accepts raw data calls, `r4vn_surv`, `r4vn_tabmulti`, `r4vn_result`, `lm`, `glm`, and `coxph` objects.", call. = FALSE)
}

.r4vn_tf_fmt_num <- function(x, digits = 2) {
  if (!length(x) || is.na(x) || !is.finite(x)) return("")
  formatC(x, format = "f", digits = digits)
}

.r4vn_tf_fmt_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_tf_text_full <- function(lang = c("en", "vi"), text = NULL) {
  lang <- match.arg(lang)
  tx <- .r4vn_tf_text(lang, text)
  extra <- if (lang == "vi") list(
    model = "M\u00f4 h\u00ecnh",
    outcome = "K\u1ebft c\u1ee5c",
    subgroup = "Ph\u00e2n nh\u00f3m",
    interaction_p = "p t\u01b0\u01a1ng t\u00e1c",
    overall = "Chung",
    effect = "\u01af\u1edbc t\u00ednh"
  ) else list(
    model = "Model",
    outcome = "Outcome",
    subgroup = "Subgroup",
    interaction_p = "P for interaction",
    overall = "Overall",
    effect = "Estimate"
  )
  for (nm in names(extra)) if (is.null(tx[[nm]])) tx[[nm]] <- extra[[nm]]
  tx
}

.r4vn_tf_publication_data <- function(rows, estimates, model_keys, display_models,
                                      effect_title, tx, effect_digit = 2,
                                      p_digit = 3, pvalue = TRUE,
                                      global_p = FALSE, show_n = FALSE,
                                      show_events = FALSE) {
  out <- data.frame(.r4vn_row_id = seq_len(nrow(rows)), stringsAsFactors = FALSE)
  out[[tx$characteristic]] <- as.character(rows$label)

  estimate_text <- function(z, reference = FALSE) {
    if (reference) return(tx$reference)
    if (!nrow(z) || !is.finite(z$estimate[1L])) return("")
    a <- paste0(.r4vn_tf_fmt_num(z$estimate[1L], effect_digit), " (",
                .r4vn_tf_fmt_num(z$lower[1L], effect_digit), "\u2013",
                .r4vn_tf_fmt_num(z$upper[1L], effect_digit), ")")
    extra <- character()
    if (isTRUE(show_n) && is.finite(z$n[1L])) extra <- c(extra, paste0(tx$n, "=", as.integer(z$n[1L])))
    if (isTRUE(show_events) && is.finite(z$events[1L])) extra <- c(extra, paste0(tx$events, "=", as.integer(z$events[1L])))
    if (length(extra)) a <- paste0(a, "; ", paste(extra, collapse = ", "))
    a
  }

  for (j in seq_along(model_keys)) {
    key <- model_keys[j]
    model_label <- as.character(display_models[[key]])
    est_col <- paste(model_label, effect_title)
    p_col <- paste(model_label, tx$p)
    est_values <- character(nrow(rows))
    p_values <- character(nrow(rows))

    for (i in seq_len(nrow(rows))) {
      r <- rows[i, , drop = FALSE]
      if (identical(r$row_type, "header")) {
        z <- estimates[estimates$variable == r$variable & estimates$model == key, , drop = FALSE]
        if (isTRUE(global_p) && nrow(z)) {
          gp <- z$global_p[is.finite(z$global_p)]
          if (length(gp)) p_values[i] <- .r4vn_tf_fmt_p(gp[1L], p_digit)
        }
      } else {
        if (identical(r$row_type, "level")) {
          z <- estimates[estimates$variable == r$variable & estimates$level == r$level & estimates$model == key, , drop = FALSE]
        } else {
          z <- estimates[estimates$variable == r$variable & (is.na(estimates$level) | estimates$level == "") & estimates$model == key, , drop = FALSE]
        }
        if (nrow(z)) {
          ref <- isTRUE(z$reference[1L])
          est_values[i] <- estimate_text(z, ref)
          if (isTRUE(pvalue) && !ref) p_values[i] <- .r4vn_tf_fmt_p(z$p[1L], p_digit)
        }
      }
    }

    out[[est_col]] <- est_values
    if (isTRUE(pvalue) || isTRUE(global_p)) out[[p_col]] <- p_values
  }

  out$.r4vn_row_id <- NULL
  names(out) <- make.unique(names(out), sep = "_")
  out
}

.r4vn_tf_resolve_named <- function(x, keys, default) {
  if (!length(keys)) return(x)
  if (length(default) == 1L) default <- rep(default, length(keys))
  default <- rep(default, length.out = length(keys))
  names(default) <- keys
  if (is.null(x)) return(default)
  if (!is.null(names(x))) {
    out <- default
    for (k in keys) if (k %in% names(x)) out[[k]] <- x[[k]]
    return(out)
  }
  out <- rep(x, length.out = length(keys))
  names(out) <- keys
  out
}

.r4vn_tf_wrap <- function(x, width) {
  if (!is.finite(width)) return(as.character(x))
  z <- strwrap(as.character(x), width = max(8L, as.integer(width)))
  if (!length(z)) "" else paste(z, collapse = "\n")
}

.r4vn_tf_expand_rows <- function(rows, model_keys, display_models,
                                 row_layout = c("auto", "modelrows", "compact"),
                                 show_model_label = TRUE) {
  row_layout <- match.arg(row_layout)
  if (identical(row_layout, "auto")) row_layout <- if (length(model_keys) > 1L) "modelrows" else "compact"
  out <- rows
  out$.display_model <- NA_character_
  out$.display_model_label <- ""
  out$.item_id <- seq_len(nrow(out))
  out$.row_group <- if (nrow(out)) cumsum(c(TRUE, out$variable[-1L] != out$variable[-nrow(out)])) else integer()
  if (identical(row_layout, "compact")) return(out)

  ans <- list()
  item_id <- 0L
  group_id <- 0L
  prev_variable <- NULL
  for (i in seq_len(nrow(rows))) {
    r <- rows[i, , drop = FALSE]
    item_id <- item_id + 1L
    if (is.null(prev_variable) || !identical(as.character(r$variable[1L]), prev_variable)) group_id <- group_id + 1L
    prev_variable <- as.character(r$variable[1L])
    if (identical(r$row_type, "header")) {
      r$.display_model <- NA_character_
      r$.display_model_label <- ""
      r$.item_id <- item_id
      r$.row_group <- group_id
      ans[[length(ans) + 1L]] <- r
      next
    }
    for (k in model_keys) {
      z <- r
      z$.display_model <- k
      z$.display_model_label <- if (isTRUE(show_model_label)) as.character(display_models[[k]]) else ""
      z$.item_id <- item_id
      z$.row_group <- group_id
      ans[[length(ans) + 1L]] <- z
    }
  }
  out <- do.call(rbind, ans)
  rownames(out) <- NULL
  out
}

.r4vn_tf_row_positions <- function(rows_plot, row_spacing = 1,
                                   model_row_gap = 0.55,
                                   group_gap = 0.25) {
  n <- nrow(rows_plot)
  if (!n) return(numeric())
  rs <- suppressWarnings(as.numeric(row_spacing)[1L]); if (!is.finite(rs) || rs <= 0) rs <- 1
  mg <- suppressWarnings(as.numeric(model_row_gap)[1L]); if (!is.finite(mg) || mg <= 0) mg <- max(.35, .55 * rs)
  gg <- suppressWarnings(as.numeric(group_gap)[1L]); if (!is.finite(gg) || gg < 0) gg <- .25 * rs

  pos <- numeric(n)
  if (n > 1L) {
    for (i in 2:n) {
      prev <- rows_plot[i - 1L, , drop = FALSE]
      cur <- rows_plot[i, , drop = FALSE]
      same_item <- isTRUE(prev$.item_id == cur$.item_id)
      gap <- if (same_item) mg else rs
      if (!same_item && !isTRUE(prev$.row_group == cur$.row_group)) gap <- gap + gg
      pos[i] <- pos[i - 1L] + gap
    }
  }
  max(pos) - pos + 1.4
}

.r4vn_tf_draw_zebra <- function(rows_plot, row_y,
                                zebra_fill = c("white", "gray94"),
                                zebra_by = c("variable", "header"),
                                xleft = 0.005, xright = 0.995) {
  zebra_by <- match.arg(zebra_by)
  if (length(zebra_fill) < 2L) zebra_fill <- rep(zebra_fill, 2L)
  group_id <- if (identical(zebra_by, "header")) rows_plot$.row_group else rows_plot$.row_group
  ug <- unique(group_id)
  ug <- ug[is.finite(ug)]
  if (!length(ug)) return(invisible(NULL))
  diffs <- diff(sort(unique(row_y)))
  pad <- if (length(diffs)) min(diffs) * .45 else .4
  for (j in seq_along(ug)) {
    idx <- which(group_id == ug[j])
    if (!length(idx)) next
    graphics::rect(xleft, min(row_y[idx]) - pad, xright, max(row_y[idx]) + pad,
                   col = zebra_fill[(j - 1L) %% length(zebra_fill) + 1L], border = NA)
  }
  invisible(NULL)
}

.r4vn_tf_numeric_text <- function(z, ref, tx, effect_digit = 2, p_digit = 3,
                                  pvalue = TRUE, show_n = FALSE,
                                  show_events = FALSE, p_layout = "inline") {
  if (ref) return(list(effect = tx$reference, p = ""))
  if (!nrow(z) || !is.finite(z$estimate[1L])) return(list(effect = "", p = ""))
  txt <- paste0(.r4vn_tf_fmt_num(z$estimate[1L], effect_digit), " (",
                .r4vn_tf_fmt_num(z$lower[1L], effect_digit), "\u2013",
                .r4vn_tf_fmt_num(z$upper[1L], effect_digit), ")")
  extra <- character()
  if (isTRUE(show_n) && is.finite(z$n[1L])) extra <- c(extra, paste0(tx$n, "=", as.integer(z$n[1L])))
  if (isTRUE(show_events) && is.finite(z$events[1L])) extra <- c(extra, paste0(tx$events, "=", as.integer(z$events[1L])))
  if (length(extra)) txt <- paste0(txt, "; ", paste(extra, collapse = ", "))
  pp <- ""
  if (isTRUE(pvalue)) pp <- .r4vn_tf_fmt_p(z$p[1L], p_digit)
  if (isTRUE(pvalue) && nzchar(pp) && identical(p_layout, "inline")) txt <- paste0(txt, "; ", tx$p, "=", pp)
  list(effect = txt, p = pp)
}

.r4vn_tf_axis_range <- function(est, null, xmin = NULL, xmax = NULL,
                                ticks = NULL, use_log = FALSE, ratio = FALSE) {
  vals <- c(est$lower, est$upper, est$estimate)
  vals <- vals[is.finite(vals)]
  if (use_log) vals <- vals[vals > 0]
  if (!length(vals)) stop("There are no finite estimates to plot.", call. = FALSE)

  if (ratio && !is.null(xmax) && is.null(xmin) && length(xmax) == 1L && is.finite(xmax) && xmax > 1) xmin <- 1 / xmax
  trans <- if (use_log) log else identity
  inv <- if (use_log) exp else identity

  if (is.null(xmin) || is.null(xmax)) {
    tv <- trans(vals)
    null_t <- trans(null)
    span <- range(c(tv, null_t), finite = TRUE)
    d <- diff(span)
    pad <- max(d * .10, if (use_log) .15 else .10)
    if (!is.finite(pad) || pad <= 0) pad <- if (use_log) .15 else .10
    if (is.null(xmin)) xmin <- inv(span[1L] - pad)
    if (is.null(xmax)) xmax <- inv(span[2L] + pad)
  }
  xmin <- as.numeric(xmin)[1L]; xmax <- as.numeric(xmax)[1L]
  if (!is.finite(xmin) || !is.finite(xmax) || xmin >= xmax) stop("`xmin` and `xmax` must define an increasing finite range.", call. = FALSE)
  if (use_log && xmin <= 0) stop("A logarithmic forest axis requires `xmin > 0`.", call. = FALSE)

  tmin <- trans(xmin); tmax <- trans(xmax)
  if (is.null(ticks)) {
    if (use_log && ratio) {
      candidate <- c(.01, .02, .05, .1, .2, .25, .5, 1, 2, 4, 5, 10, 20, 50, 100)
      ticks <- candidate[candidate >= xmin & candidate <= xmax]
      if (length(ticks) > 7L) {
        keep <- unique(round(seq(1, length(ticks), length.out = 7L)))
        ticks <- ticks[keep]
        if (xmin <= 1 && xmax >= 1 && !1 %in% ticks) ticks <- sort(unique(c(ticks, 1)))
      }
      if (length(ticks) < 2L) {
        ticks <- inv(pretty(c(tmin, tmax), n = 5))
      }
    } else if (use_log) {
      ticks <- inv(pretty(c(tmin, tmax), n = 5))
    } else {
      ticks <- pretty(c(xmin, xmax), n = 5)
    }
  }
  ticks <- as.numeric(ticks)
  ticks <- ticks[is.finite(ticks) & ticks >= xmin & ticks <= xmax & (!use_log | ticks > 0)]
  list(xmin = xmin, xmax = xmax, ticks = ticks, trans = trans, inv = inv,
       tmin = tmin, tmax = tmax)
}

.r4vn_tf_device <- function(file, width, height, dpi) {
  ext <- tolower(tools::file_ext(file))
  if (ext == "pdf") {
    if (isTRUE(capabilities("cairo"))) {
      grDevices::cairo_pdf(file, width = width, height = height, onefile = TRUE)
    } else {
      grDevices::pdf(file, width = width, height = height, onefile = TRUE)
    }
  }
  else if (ext == "png") grDevices::png(file, width = width, height = height, units = "in", res = dpi)
  else if (ext == "svg") grDevices::svg(file, width = width, height = height)
  else if (ext %in% c("jpg", "jpeg")) grDevices::jpeg(file, width = width, height = height, units = "in", res = dpi, quality = 95)
  else if (ext %in% c("tif", "tiff")) grDevices::tiff(file, width = width, height = height, units = "in", res = dpi, compression = "lzw")
  else stop("Unsupported graphics extension. Use pdf, png, svg, jpg/jpeg, or tiff.", call. = FALSE)
}

.r4vn_tf_panel_value <- function(x, panel_name, panel_index, default = NULL) {
  if (is.null(x)) return(default)
  if (is.list(x)) {
    if (!is.null(names(x)) && panel_name %in% names(x)) return(x[[panel_name]])
    if (length(x) >= panel_index) return(x[[panel_index]])
    return(default)
  }
  if (!is.null(names(x)) && panel_name %in% names(x)) return(x[[panel_name]])
  if (length(x) == 1L) return(x[[1L]])
  if (length(x) >= panel_index) return(x[[panel_index]])
  default
}

.r4vn_tf_outcome_specs <- function(outcomes, data) {
  if (is.null(outcomes)) return(NULL)
  if (is.character(outcomes)) {
    vals <- as.character(outcomes)
    labs <- names(outcomes)
    if (is.null(labs) || any(!nzchar(labs))) labs <- vals
    out <- lapply(seq_along(vals), function(i) list(outcome = vals[i], label = labs[i]))
    names(out) <- labs
    return(out)
  }
  if (!is.list(outcomes) || !length(outcomes)) stop("`outcomes` must be a named character vector or a named list of outcome specifications.", call. = FALSE)
  labs <- names(outcomes)
  if (is.null(labs)) labs <- rep("", length(outcomes))
  out <- vector("list", length(outcomes))
  for (i in seq_along(outcomes)) {
    z <- outcomes[[i]]
    if (is.character(z) && length(z) == 1L) z <- list(outcome = z)
    if (!is.list(z) || is.null(z$outcome)) stop("Each element of `outcomes` must contain `outcome`.", call. = FALSE)
    z$outcome <- as.character(z$outcome)[1L]
    if (!z$outcome %in% names(data)) stop(sprintf("Outcome `%s` was not found in `data`.", z$outcome), call. = FALSE)
    if (!is.null(z$time)) {
      z$time <- as.character(z$time)[1L]
      if (!z$time %in% names(data)) stop(sprintf("Time variable `%s` was not found in `data`.", z$time), call. = FALSE)
    }
    lab <- if (nzchar(labs[i])) labs[i] else if (!is.null(z$label)) as.character(z$label)[1L] else z$outcome
    z$label <- lab
    out[[i]] <- z
  }
  names(out) <- vapply(out, function(z) z$label, character(1))
  out
}

.r4vn_tf_formula_interaction <- function(response, predictor, subgroup, covariates = character()) {
  p <- .r4vn_tf_escape_name(predictor)
  s <- .r4vn_tf_escape_name(subgroup)
  rhs <- paste0(p, " * ", s)
  if (length(covariates)) rhs <- paste(rhs, paste(vapply(covariates, .r4vn_tf_escape_name, character(1)), collapse = " + "), sep = " + ")
  stats::as.formula(paste(response, "~", rhs))
}

.r4vn_tf_fit_interaction <- function(d, predictor, subgroup, covariates, effect) {
  response <- if (identical(effect, "HR")) "survival::Surv(.time, .status)" else ".outcome"
  f <- .r4vn_tf_formula_interaction(response, predictor, subgroup, covariates)
  if (identical(effect, "HR")) {
    if (!requireNamespace("survival", quietly = TRUE)) stop("Cox subgroup forests require the `survival` package.", call. = FALSE)
    fit <- survival::coxph(f, data = d, x = TRUE, y = TRUE, model = TRUE, ties = "efron")
    V <- stats::vcov(fit)
  } else if (identical(effect, "Beta")) {
    fit <- stats::lm(f, data = d, x = TRUE, y = TRUE)
    V <- stats::vcov(fit)
  } else if (identical(effect, "OR")) {
    fit <- stats::glm(f, family = stats::binomial("logit"), data = d, x = TRUE, y = TRUE)
    V <- stats::vcov(fit)
  } else {
    fit <- stats::glm(f, family = stats::poisson("log"), data = d, x = TRUE, y = TRUE)
    V <- if (effect %in% c("RR", "PR")) .r4vn_tf_robust_poisson(fit) else stats::vcov(fit)
  }
  list(fit = fit, vcov = V)
}

.r4vn_tf_interaction_p <- function(fit_obj, predictor, subgroup) {
  fit <- fit_obj$fit; V <- fit_obj$vcov
  mm <- stats::model.matrix(fit)
  ass <- attr(mm, "assign")
  tl <- attr(stats::terms(fit), "term.labels")
  clean <- gsub("`", "", tl, fixed = TRUE)
  int_idx <- which(vapply(clean, function(z) {
    parts <- strsplit(z, ":", fixed = TRUE)[[1L]]
    predictor %in% parts && subgroup %in% parts
  }, logical(1)))
  if (!length(int_idx)) return(NA_real_)
  cn <- colnames(mm)[ass %in% int_idx]
  b <- stats::coef(fit)
  cn <- intersect(cn, names(b))
  cn <- cn[is.finite(b[cn])]
  if (!length(cn)) return(NA_real_)
  vv <- V[cn, cn, drop = FALSE]
  good <- is.finite(diag(vv)) & diag(vv) > 0
  cn <- cn[good]
  if (!length(cn)) return(NA_real_)
  bb <- b[cn]; vv <- V[cn, cn, drop = FALSE]
  inv <- tryCatch(solve(vv), error = function(e) tryCatch(qr.solve(vv), error = function(e2) NULL))
  if (is.null(inv)) return(NA_real_)
  q <- as.numeric(t(bb) %*% inv %*% bb)
  df <- qr(vv)$rank
  if (!is.finite(q) || df < 1L) return(NA_real_)
  stats::pchisq(q, df = df, lower.tail = FALSE)
}

.r4vn_tf_subgroup_adjust_spec <- function(data, predictor_spec, adjusted, multi) {
  if (isTRUE(adjusted)) stop("In subgroup mode, use `adjusted = vars(...)` rather than `adjusted = TRUE` so the adjustment set is explicit.", call. = FALSE)
  out <- NULL
  if (!isFALSE(adjusted) && !is.null(adjusted)) out <- .r4vn_tf_spec(adjusted, data, "adjusted")
  if (!isFALSE(multi) && !is.null(multi) && !isTRUE(multi)) out <- .r4vn_tf_spec(multi, data, "multi")
  if (!is.null(out) && nrow(out)) out <- out[out$variable != predictor_spec$variable[1L], , drop = FALSE]
  out
}

.r4vn_tf_subgroup_build <- function(data, outcome_name, time_name,
                                    predictor_spec, subgroup_spec,
                                    adjusted = FALSE, multi = FALSE,
                                    effect, event = NULL, failure = NULL,
                                    ci = .95, sample = c("auto", "common", "model"),
                                    per = NULL, labels = NULL, level_labels = NULL) {
  sample <- match.arg(sample)
  pred <- predictor_spec$variable[1L]
  if (identical(predictor_spec$type[1L], "categorical") && length(.r4vn_tf_levels(data[[pred]])) != 2L) {
    stop("Subgroup forest currently requires a continuous predictor or a two-level categorical predictor.", call. = FALSE)
  }
  for (i in seq_len(nrow(subgroup_spec))) {
    sg <- subgroup_spec$variable[i]
    if (identical(subgroup_spec$type[i], "mean")) stop(sprintf("Subgroup variable `%s` must be categorical. Create a grouped variable first.", sg), call. = FALSE)
  }
  adj_spec <- .r4vn_tf_subgroup_adjust_spec(data, predictor_spec, adjusted, multi)
  status_name <- if (identical(effect, "HR")) outcome_name else NULL

  analysis_data <- data
  common <- identical(sample, "common")
  if (common) {
    need <- unique(c(outcome_name, time_name, pred, subgroup_spec$variable,
                     if (!is.null(adj_spec)) adj_spec$variable else NULL))
    keep <- stats::complete.cases(data[, need, drop = FALSE])
    analysis_data <- data[keep, , drop = FALSE]
  }

  est_list <- list(); rows <- list(); models <- list(); interactions <- list()
  for (i in seq_len(nrow(subgroup_spec))) {
    sg <- subgroup_spec$variable[i]
    sg_label <- .r4vn_tf_var_label(data, sg, labels)
    rows[[length(rows) + 1L]] <- data.frame(variable = sg, level = NA_character_, row_type = "header", label = sg_label, stringsAsFactors = FALSE)

    cov_spec <- adj_spec
    if (!is.null(cov_spec) && nrow(cov_spec)) cov_spec <- cov_spec[!cov_spec$variable %in% c(pred, sg), , drop = FALSE]
    full_spec <- .r4vn_tf_merge_spec(predictor_spec, subgroup_spec[i, , drop = FALSE], cov_spec)
    d_int <- tryCatch(.r4vn_tf_prepare(analysis_data, outcome_name, time_name, status_name, full_spec, effect, event, failure), error = function(e) NULL)
    ip <- NA_real_; int_fit <- NULL
    if (!is.null(d_int) && nrow(d_int)) {
      int_fit <- tryCatch(.r4vn_tf_fit_interaction(d_int, pred, sg, if (is.null(cov_spec)) character() else cov_spec$variable, effect), error = function(e) NULL)
      if (!is.null(int_fit)) ip <- .r4vn_tf_interaction_p(int_fit, pred, sg)
    }
    interactions[[sg]] <- int_fit

    levs <- .r4vn_tf_levels(data[[sg]])
    for (lv in levs) {
      rows[[length(rows) + 1L]] <- data.frame(variable = sg, level = as.character(lv), row_type = "level",
                                              label = .r4vn_tf_level_label(sg, as.character(lv), level_labels), stringsAsFactors = FALSE)
      sub <- analysis_data[as.character(analysis_data[[sg]]) == as.character(lv) & !is.na(analysis_data[[sg]]), , drop = FALSE]
      ms <- .r4vn_tf_merge_spec(predictor_spec, cov_spec)
      z <- NULL; fit <- NULL
      if (nrow(sub)) {
        z <- tryCatch({
          dd <- .r4vn_tf_prepare(sub, outcome_name, time_name, status_name, ms, effect, event, failure)
          m <- .r4vn_tf_fit(dd, ms, effect, ci)
          fit <- m$fit
          e <- .r4vn_tf_one_effect(m, pred, predictor_spec, effect, .r4vn_tf_per(pred, per), ci, "Subgroup")
          if (identical(predictor_spec$type[1L], "categorical")) e <- e[!e$reference & is.finite(e$estimate), , drop = FALSE]
          else e <- e[is.finite(e$estimate), , drop = FALSE]
          if (nrow(e)) e[1L, , drop = FALSE] else NULL
        }, error = function(e) NULL)
      }
      models[[paste(sg, lv, sep = "::")]] <- fit
      if (is.null(z)) {
        z <- data.frame(variable = sg, level = as.character(lv), reference = FALSE,
                        model = "Subgroup", estimate = NA_real_, lower = NA_real_, upper = NA_real_,
                        p = NA_real_, global_p = NA_real_, n = nrow(sub), events = NA_real_, stringsAsFactors = FALSE)
      } else {
        z$variable <- sg
        z$level <- as.character(lv)
      }
      z$interaction_p <- ip
      z$subgroup_variable <- sg
      z$subgroup_level <- as.character(lv)
      est_list[[length(est_list) + 1L]] <- z
    }
  }

  estimates <- do.call(rbind, est_list)
  rows <- do.call(rbind, rows)
  rownames(estimates) <- NULL; rownames(rows) <- NULL
  list(estimates = estimates, rows = rows, models = models, interactions = interactions,
       adjusted_spec = adj_spec, common_sample = common)
}

.r4vn_tf_subgroup_publication <- function(rows, estimates, tx, effect_title,
                                          pvalue = TRUE, show_interaction_p = TRUE,
                                          effect_digit = 2, p_digit = 3,
                                          show_n = FALSE, show_events = FALSE) {
  out <- data.frame(.id = seq_len(nrow(rows)), stringsAsFactors = FALSE)
  out[[tx$subgroup]] <- rows$label
  out[[effect_title]] <- ""
  if (isTRUE(pvalue)) out[[tx$p]] <- ""
  if (isTRUE(show_interaction_p)) out[[tx$interaction_p]] <- ""
  for (i in seq_len(nrow(rows))) {
    r <- rows[i, , drop = FALSE]
    if (identical(r$row_type, "header")) {
      if (isTRUE(show_interaction_p)) {
        z <- estimates[estimates$subgroup_variable == r$variable, , drop = FALSE]
        ip <- z$interaction_p[is.finite(z$interaction_p)]
        if (length(ip)) out[[tx$interaction_p]][i] <- .r4vn_tf_fmt_p(ip[1L], p_digit)
      }
      next
    }
    z <- estimates[estimates$subgroup_variable == r$variable & estimates$subgroup_level == r$level, , drop = FALSE]
    if (!nrow(z) || !is.finite(z$estimate[1L])) next
    txt <- paste0(.r4vn_tf_fmt_num(z$estimate[1L], effect_digit), " (",
                  .r4vn_tf_fmt_num(z$lower[1L], effect_digit), "\u2013",
                  .r4vn_tf_fmt_num(z$upper[1L], effect_digit), ")")
    extra <- character()
    if (isTRUE(show_n) && is.finite(z$n[1L])) extra <- c(extra, paste0(tx$n, "=", as.integer(z$n[1L])))
    if (isTRUE(show_events) && is.finite(z$events[1L])) extra <- c(extra, paste0(tx$events, "=", as.integer(z$events[1L])))
    if (length(extra)) txt <- paste0(txt, "; ", paste(extra, collapse = ", "))
    out[[effect_title]][i] <- txt
    if (isTRUE(pvalue)) out[[tx$p]][i] <- .r4vn_tf_fmt_p(z$p[1L], p_digit)
  }
  out$.id <- NULL
  out
}


#' Flexible regression, multi-outcome, survival, and subgroup forest plots
#'
#' `tabforest()` is the common forest-plot engine for R4VN. It can (1) fit
#' regression models directly from an outcome and focal predictors, (2) reuse a
#' fitted R4VN or standard R model, (3) place several outcomes side by side using
#' the same predictor structure, and (4) create subgroup-effect forests with a
#' p-value for interaction. Numeric estimates are always stored without clipping;
#' `xmin` and `xmax` affect only the drawing.
#'
#' @param outcome Outcome variable for ordinary regression/subgroup analysis, or a
#'   supported fitted object (`r4vn_surv`, `r4vn_tabmulti`, `r4vn_stat`, `lm`,
#'   `glm`, `coxph`). For Cox regression, `outcome` is the event/status variable.
#'   May be omitted when `outcomes` is supplied.
#' @param predictors Focal predictors that should appear in the forest. Prefer
#'   `vars()` so R4VN type/reference declarations are retained, for example
#'   `vars(c.age, b2.sex, c.bmi, smoking)`.
#' @param data Data frame. If omitted, the active R4VN data frame is used.
#' @param time Follow-up time variable for Cox regression. Supplying `time`
#'   automatically selects HR unless another incompatible effect is requested.
#' @param event Modeled event level for binary OR/RR/PR outcomes.
#' @param failure Event value for Cox regression. With numeric 0/1 status, 1 is
#'   selected automatically when present.
#' @param outcomes Optional named character vector or named list for a
#'   multi-outcome forest. Character example: `c(HTN="hypertension",
#'   DEP="depression")`. A list allows outcome-specific settings, e.g.
#'   `list(HTN=list(outcome="hypertension",event="Yes"),
#'   Death=list(outcome="death",time="followup",failure=1,estimate="hr"))`.
#' @param subgroup Optional categorical variables for subgroup analysis. When
#'   supplied, use `predictor` for the main exposure whose effect is estimated
#'   within each subgroup level.
#' @param predictor Main exposure for subgroup mode. It may be continuous or a
#'   two-level categorical variable. Use `vars(b2.treatment)` when a non-default
#'   reference is required.
#' @param type Analysis mode: `"auto"`, `"regression"`, `"multioutcome"`, or
#'   `"subgroup"`. `"auto"` chooses multi-outcome when `outcomes` is non-NULL,
#'   subgroup when `subgroup` is non-NULL, otherwise ordinary regression.
#' @param crude For regression/multi-outcome mode, fit one crude model per focal
#'   predictor. Set `FALSE` when only adjusted/multivariable estimates are wanted.
#' @param adjusted Regression mode: `FALSE`, `TRUE`, or `vars(...)`. With
#'   `vars(X1,X2)`, each focal predictor gets a separate model adjusted for X1/X2.
#'   With `TRUE`, each focal predictor is adjusted for all other focal predictors.
#'   Subgroup mode requires an explicit `vars(...)` adjustment set if adjustment
#'   is desired.
#' @param multi Regression mode: `FALSE`, `TRUE`, or `vars(...)`. `TRUE` fits one
#'   joint model containing all focal predictors. `vars(A,B,C,X)` fits that exact
#'   joint model but still displays only variables listed in `predictors`.
#'   In subgroup mode, an explicit `vars(...)` may also be used as the final
#'   adjustment set; the exposure and current subgroup variable are removed from
#'   the covariate set automatically.
#' @param or,rr,pr,irr Logical shortcuts for OR, RR, PR, or IRR. Only one may be
#'   TRUE. Binary outcomes default to OR; numeric outcomes default to beta.
#' @param estimate Explicit effect type: `"auto"`, `"beta"`, `"or"`, `"rr"`,
#'   `"pr"`, `"irr"`, or `"hr"`.
#' @param ci Confidence level, default 0.95.
#' @param sample Missing-data strategy. `"common"` forces displayed regression
#'   models to use the same complete-case sample; `"model"` allows each model to
#'   use its own available cases; `"auto"` uses a common sample when several
#'   regression model groups are displayed. In subgroup mode, `"auto"` behaves
#'   like model-specific analysis so unrelated subgroup variables do not reduce
#'   one another's sample size.
#' @param per Optional multiplier for continuous effects. Example `c(age=10)`
#'   reports the ratio/HR per 10 years or beta per 10 units.
#' @param per_labels Optional display labels for `per`, e.g.
#'   `c(age="per 10 years")`.
#' @param select Model components when `outcome` is a fitted R4VN object. For
#'   `tabsurv()` this can include `"crude"`, `"adjusted"`, `"multi"`; for
#'   `tabmulti()` use stored model names such as `"full"`, `"backward"`.
#' @param xmin,xmax Forest plotting limits. These never alter stored estimates.
#'   For ratio effects, if only `xmax` is supplied and is >1, `xmin=1/xmax` is
#'   used automatically. In multi-outcome mode these may be named vectors or
#'   lists keyed by outcome-panel name.
#' @param ticks Optional axis ticks. In multi-outcome mode a named list can give
#'   different ticks to different outcome panels.
#' @param log `NULL` or logical. Ratio effects default to logarithmic axes, beta
#'   to linear. In multi-outcome mode this may be a named logical vector/list.
#' @param arrows Draw arrowheads when CIs extend beyond plotting limits. When the
#'   point estimate itself is outside the range, no false boundary point is drawn.
#' @param row_layout `"auto"`, `"modelrows"`, or `"compact"`. The default
#'   `"auto"` uses `"modelrows"` whenever more than one model group is displayed
#'   and `"compact"` when only one model is displayed. In `"modelrows"`, Crude,
#'   Adjusted, and Multivariable are separate physical rows but share ONE effect
#'   column (for example, one `OR (95% CI)` column). Thus each CI, marker, numeric
#'   estimate, and p-value is aligned with its own row. Use `"compact"` only when
#'   several model estimates are deliberately wanted on the same labelled row.
#' @param layout In compact mode, `"dodge"` or `"stack"` controls vertical
#'   offsets of multiple model markers/CI lines.
#' @param row_spacing Baseline distance between ordinary rows.
#' @param model_row_gap Distance between model rows belonging to the same
#'   variable/level when `row_layout="modelrows"`.
#' @param group_gap Extra vertical separation between variable blocks.
#' @param order Optional order of focal predictor variable names.
#' @param reference Show categorical reference rows. Default TRUE.
#' @param pvalue Show coefficient-level p-values.
#' @param global_p Show categorical-variable omnibus Wald p-values. Default
#'   FALSE because forest plots are usually cleaner without these values.
#' @param show_n,show_events Add model N and number of events after the estimate.
#' @param show_model_label In `modelrows`, print the model label (Crude, Adjusted,
#'   Multivariable) beside the corresponding row.
#' @param show_interaction_p In subgroup mode, show the p-value for interaction on
#'   the subgroup-variable header row.
#' @param p_layout Regression display style: `"inline"` appends p to the numeric
#'   effect string; `"column"` uses a separate p-value column. `modelrows` works
#'   particularly well with either style.
#' @param effect_digit,p_digit Decimal places for effects and p-values.
#' @param labels Named character vector overriding variable labels.
#' @param level_labels Named list overriding displayed categorical levels. This
#'   changes display only, not model coding/reference levels.
#' @param lang Built-in language: `"en"` or `"vi"`.
#' @param text Named list overriding individual words. Useful keys include
#'   `characteristic`, `reference`, `crude`, `adjusted`, `multi`, `p`, `n`,
#'   `events`, `subgroup`, `interaction_p`, `overall`, `arrow_note`, and effect
#'   keys `OR`, `RR`, `PR`, `IRR`, `HR`, `beta`.
#' @param model_labels Named character vector overriding model labels, e.g.
#'   `c(Crude="Unadjusted",Multivariable="Adjusted")`.
#' @param label_title,effect_title,axis_title Optional column/axis titles.
#' @param title,subtitle,caption Optional plot title, subtitle, caption.
#' @param note TRUE for an automatic clipping note, FALSE for none, or custom text.
#' @param template Visual preset: `"journal"`, `"clean"`, `"minimal"`.
#' @param grid `"major"`, `"none"`, or `"both"`.
#' @param font_family Base graphics font family.
#' @param base_size,label_cex,header_cex,axis_cex,model_cex Text-size controls.
#' @param colors Model line/marker colors. A named vector is recommended, e.g.
#'   `c(Crude="gray50",Adjusted="navy",Multivariable="firebrick")`.
#' @param fills Optional marker fill colors, useful with pch 21:25.
#' @param pch Model marker symbols. Named vectors may assign different symbols to
#'   crude and adjusted estimates.
#' @param lty Model CI line types.
#' @param point_cex Model marker sizes. May be scalar or named vector by model.
#' @param point_lwd Marker border widths. May be scalar or named vector.
#' @param ci_lwd CI line widths. May be scalar or named vector by model.
#' @param ref_lwd,ref_lty,ref_col Null-line appearance.
#' @param arrow_length Arrowhead size in inches.
#' @param zebra Draw alternating background blocks by predictor/subgroup.
#' @param zebra_fill Two or more background colors, e.g.
#'   `c("white","gray93")`.
#' @param zebra_by `"variable"` or `"header"`; both currently alternate complete
#'   variable/subgroup blocks so all levels/models in a block share a background.
#' @param label_width,forest_width,column_gap Horizontal layout controls for a
#'   single regression/subgroup forest.
#' @param panel_gap Gap between panels in multi-outcome mode.
#' @param panel_forest_ratio Fraction of each multi-outcome panel devoted to the
#'   CI forest; the remaining panel width is used for numeric estimates.
#' @param label_indent Indentation of categorical levels.
#' @param label_wrap Approximate wrapping width for long labels; Inf disables.
#' @param file Optional PDF/PNG/SVG/JPG/TIFF output file.
#' @param width,height,dpi Graphics dimensions. If height is NULL it grows with
#'   the actual number of drawn rows, so `modelrows` can produce a tall figure
#'   without compressing row spacing.
#' @param show Draw immediately. Default TRUE.
#' @param console Print the long standardized estimate table.
#'
#' @return An object of class `r4vn_tabforest`; multi-outcome and subgroup modes
#'   add subclasses `r4vn_tabforest_multi` and `r4vn_tabforest_subgroup`.
#'   `$data` is publication-ready, `$table` is the numeric long table, `$models`
#'   stores fitted models, and `plot()` can redraw without refitting.
#'
#' @details
#' ## 1. Regression model semantics
#'
#' With `predictors=vars(A,B,C)`:
#'
#' * `crude=TRUE`: Y~A, Y~B, Y~C.
#' * `adjusted=vars(X)`: Y~A+X, Y~B+X, Y~C+X.
#' * `adjusted=TRUE`: each focal variable is adjusted for the other focal vars.
#' * `multi=TRUE`: one joint model Y~A+B+C.
#' * `multi=vars(A,B,C,X)`: one joint model Y~A+B+C+X, but only A/B/C are shown.
#'
#' Therefore `adjusted` and `multi` answer different scientific questions and may
#' be requested together in the same forest.
#' When two or more model groups are displayed, `row_layout="auto"` uses separate
#' model rows and ONE shared effect column. For example, Crude and Multivariable
#' ORs are both printed under the same `OR (95% CI)` header instead of being put
#' in separate Crude-OR and Multivariable-OR columns.
#'
#' ## 2. Effect measure selected by outcome
#'
#' * numeric continuous outcome -> linear regression -> beta, null=0;
#' * binary outcome -> logistic regression -> OR, null=1;
#' * `pr=TRUE` -> robust modified Poisson -> PR, null=1;
#' * `rr=TRUE` -> robust modified Poisson -> RR, null=1;
#' * count outcome + `irr=TRUE` -> Poisson -> IRR, null=1;
#' * `time=` -> Cox proportional hazards -> HR, null=1.
#'
#' ## 3. Exact numbers when the forest is clipped
#'
#' Suppose OR=7.41 and 95% CI=1.56 to 35.20 while `xmax=10`. The printed number
#' remains `7.41 (1.56-35.20)`. Only the graphical CI is truncated at 10 and an
#' arrow is drawn. If OR itself exceeds 10, no marker is placed falsely at 10.
#'
#' ## 4. Multiple models: separate rows, one merged effect column
#'
#' `row_layout="auto"` is the default. If more than one model group is present,
#' it automatically switches to the `"modelrows"` layout. Crude, Adjusted and/or
#' Multivariable estimates are placed on separate physical rows, while the right
#' side contains only ONE shared effect column such as `OR (95% CI)`, `HR (95% CI)`,
#' or `Beta (95% CI)`. The result, CI line, marker and p-value therefore stay on
#' exactly the same row. This is the recommended publication layout when crude and
#' adjusted estimates are presented together. Increase `row_spacing`,
#' `model_row_gap`, `group_gap`, or leave `height=NULL` for a taller figure.
#'
#' Set `row_layout="compact"` only when you intentionally want several model
#' estimates on the same labelled row; compact mode retains separate numeric model
#' columns because the rows are not expanded.
#'
#' ## 5. Multi-outcome forests
#'
#' `outcomes=` creates side-by-side panels sharing predictor labels. Every panel
#' may have its own effect type, axis range and follow-up variable. This permits
#' two binary outcomes (OR panels), several continuous outcomes (beta panels), or
#' even mixed OR/HR panels in one figure. For very many panels, increase `width`.
#'
#' ## 6. Subgroup forests
#'
#' `subgroup=` estimates the effect of one main `predictor` separately within
#' each subgroup level. A full model containing predictor*subgroup is fitted for
#' the Wald interaction p-value. For RR/PR, the interaction Wald test uses the
#' same robust covariance approach as the effect model. Cox subgroup forests use
#' HR and a Cox interaction model. A subgroup variable must be categorical; make
#' clinically meaningful groups before calling `tabforest()`.
#'
#' ## 7. Styling
#'
#' `colors`, `fills`, `pch`, `lty`, `point_cex`, `point_lwd`, and `ci_lwd` accept
#' named model vectors. `zebra=TRUE` shades complete variable blocks, closely
#' matching journal forest-table layouts. `lang="vi"`, `text=`, `labels=`, and
#' `level_labels=` allow all visible wording to be translated without changing
#' the model.
#'
#' @seealso `vars`, `tab`, `tabmulti`, `tabsurv`, `tabexport`
#' @family R4VN tables
#'
#' @examples
#' data(tabforest_demo)
#' usedf(tabforest_demo, quiet = TRUE)
#'
#' # Crude odds ratios.
#' f1 <- tabforest(
#'   hypertension,
#'   predictors = vars(c.age, sex, c.bmi, smoking),
#'   event = "Yes",
#'   show = FALSE
#' )
#' f1$data
#'
#' # Crude plus one final multivariable model.
#' f2 <- tabforest(
#'   hypertension,
#'   predictors = vars(c.age, sex, c.bmi, smoking),
#'   event = "Yes",
#'   crude = TRUE, multi = TRUE,
#'   show = FALSE
#' )
#'
#' # Re-drawing is intentionally interactive so CRAN examples do not depend
#' # on the graphics device or installed fonts.
#' if (interactive()) {
#'   plot(f2, row_layout = "modelrows", zebra = TRUE)
#' }
#'
#' # Each focal predictor adjusted for the same confounders.
#' f3 <- tabforest(
#'   hypertension,
#'   predictors = vars(c.age, c.bmi, smoking),
#'   event = "Yes",
#'   adjusted = vars(sex, education),
#'   show = FALSE
#' )
#'
#' # Vietnamese display text can be prepared without drawing during checks.
#' f_vi <- tabforest(
#'   hypertension,
#'   predictors = vars(c.age, sex, c.bmi, smoking),
#'   event = "Yes",
#'   multi = TRUE,
#'   lang = "vi",
#'   labels = c(
#'     age = "Tu\u1ed5i",
#'     sex = "Gi\u1edbi t\u00ednh",
#'     bmi = "Ch\u1ec9 s\u1ed1 kh\u1ed1i c\u01a1 th\u1ec3",
#'     smoking = "H\u00fat thu\u1ed1c"
#'   ),
#'   level_labels = list(
#'     sex = c(Female = "N\u1eef", Male = "Nam"),
#'     smoking = c(No = "Kh\u00f4ng", Yes = "C\u00f3")
#'   ),
#'   text = list(reference = "Tham chi\u1ebfu"),
#'   title = "Bi\u1ec3u \u0111\u1ed3 forest",
#'   show = FALSE
#' )
#' if (interactive()) plot(f_vi)
#'
#' \donttest{
#' # Modified-Poisson prevalence ratio.
#' f_pr <- tabforest(
#'   depression,
#'   predictors = vars(c.age, sex, smoking, alcohol),
#'   event = "Yes", pr = TRUE, multi = TRUE,
#'   show = FALSE
#' )
#'
#' # Continuous outcome.
#' f_beta <- tabforest(
#'   sbp,
#'   predictors = vars(c.age, sex, c.bmi, smoking),
#'   multi = TRUE,
#'   show = FALSE
#' )
#'
#' # Cox model, only when the suggested package is available.
#' if (requireNamespace("survival", quietly = TRUE)) {
#'   f_hr <- tabforest(
#'     death, time = followup,
#'     predictors = vars(c.age, sex, treatment, c.bmi),
#'     failure = 1, crude = TRUE, multi = TRUE,
#'     show = FALSE
#'   )
#' }
#'
#' # Multi-outcome forest without drawing.
#' f_multi <- tabforest(
#'   outcomes = c(
#'     Hypertension = "hypertension",
#'     Depression = "depression"
#'   ),
#'   predictors = vars(c.age, sex, c.bmi, smoking),
#'   event = "Yes", crude = FALSE, multi = TRUE,
#'   show = FALSE
#' )
#'
#' # Subgroup forest without drawing.
#' f_sub <- tabforest(
#'   hypertension,
#'   predictor = vars(treatment),
#'   subgroup = vars(age_group, sex, obesity, diabetes, smoking),
#'   event = "Yes", type = "subgroup",
#'   adjusted = vars(c.age, c.bmi),
#'   show = FALSE
#' )
#'
#' # File output uses a temporary path and is cleaned up.
#' f_png <- tempfile(fileext = ".png")
#' tabforest(
#'   hypertension,
#'   predictors = vars(c.age, sex, c.bmi, smoking),
#'   event = "Yes", multi = TRUE,
#'   file = f_png, width = 8, height = 5, dpi = 120,
#'   show = FALSE
#' )
#' unlink(f_png)
#' }
#'
#' usedf(clear = TRUE, quiet = TRUE)
#' @export
tabforest <- function(outcome = NULL, predictors = NULL, data = NULL,
                      time = NULL, event = NULL, failure = NULL,
                      outcomes = NULL, subgroup = NULL, predictor = NULL,
                      type = c("auto", "regression", "multioutcome", "subgroup"),
                      crude = TRUE, adjusted = FALSE, multi = FALSE,
                      or = FALSE, rr = FALSE, pr = FALSE, irr = FALSE,
                      estimate = c("auto", "beta", "or", "rr", "pr", "irr", "hr"),
                      ci = 0.95, sample = c("auto", "common", "model"),
                      per = NULL, per_labels = NULL, select = NULL,
                      xmin = NULL, xmax = NULL, ticks = NULL, log = NULL,
                      arrows = TRUE,
                      row_layout = c("auto", "modelrows", "compact"),
                      layout = c("dodge", "stack"),
                      row_spacing = 1, model_row_gap = 0.55, group_gap = 0.25,
                      order = NULL, reference = TRUE,
                      pvalue = TRUE, global_p = FALSE,
                      show_n = FALSE, show_events = FALSE,
                      show_model_label = TRUE, show_interaction_p = TRUE,
                      p_layout = c("inline", "column"),
                      effect_digit = 2, p_digit = 3,
                      labels = NULL, level_labels = NULL,
                      lang = c("en", "vi"), text = NULL, model_labels = NULL,
                      label_title = NULL, effect_title = NULL, axis_title = NULL,
                      title = NULL, subtitle = NULL, caption = NULL, note = TRUE,
                      template = c("journal", "clean", "minimal"),
                      grid = c("major", "none", "both"), font_family = "",
                      base_size = 11, label_cex = 1, header_cex = 1,
                      axis_cex = 1, model_cex = 0.86,
                      colors = NULL, fills = NULL, pch = NULL, lty = NULL,
                      point_cex = 1.15, point_lwd = 1, ci_lwd = 1.2,
                      ref_lwd = 1, ref_lty = 2, ref_col = "gray45",
                      arrow_length = 0.08,
                      zebra = FALSE, zebra_fill = c("white", "gray94"),
                      zebra_by = c("variable", "header"),
                      label_width = 0.28, forest_width = 0.32,
                      column_gap = 0.012, panel_gap = 0.012,
                      panel_forest_ratio = 0.58,
                      label_indent = 0.018, label_wrap = 38,
                      file = NULL, width = 12, height = NULL, dpi = 300,
                      show = TRUE, console = FALSE) {
  call <- match.call()
  env <- parent.frame()
  type <- match.arg(type)
  lang <- match.arg(lang)
  row_layout <- match.arg(row_layout)
  layout <- match.arg(layout)
  template <- match.arg(template)
  grid <- match.arg(grid)
  zebra_by <- match.arg(zebra_by)
  p_layout <- match.arg(p_layout)
  sample <- match.arg(sample)
  if (type == "auto") type <- if (!is.null(outcomes)) "multioutcome" else if (!is.null(subgroup)) "subgroup" else "regression"
  if (!is.numeric(ci) || length(ci) != 1L || !is.finite(ci) || ci <= 0 || ci >= 1) stop("`ci` must be between 0 and 1.", call. = FALSE)
  for (nm in c("effect_digit", "p_digit")) {
    val <- get(nm)
    if (!is.numeric(val) || length(val) != 1L || is.na(val) || val < 0 || val != floor(val)) stop(sprintf("`%s` must be a non-negative integer.", nm), call. = FALSE)
  }

  common_settings <- list(
    ci=ci, xmin=xmin, xmax=xmax, ticks=ticks, log=log, arrows=arrows,
    row_layout=row_layout, layout=layout, row_spacing=row_spacing,
    model_row_gap=model_row_gap, group_gap=group_gap,
    reference=reference, pvalue=pvalue, global_p=global_p,
    show_n=show_n, show_events=show_events,
    show_model_label=show_model_label, show_interaction_p=show_interaction_p,
    p_layout=p_layout, effect_digit=effect_digit, p_digit=p_digit,
    labels=labels, level_labels=level_labels, lang=lang, text=text,
    model_labels=model_labels, label_title=label_title,
    effect_title=effect_title, axis_title=axis_title,
    title=title, subtitle=subtitle, caption=caption, note=note,
    template=template, grid=grid, font_family=font_family,
    base_size=base_size, label_cex=label_cex, header_cex=header_cex,
    axis_cex=axis_cex, model_cex=model_cex,
    colors=colors, fills=fills, pch=pch, lty=lty,
    point_cex=point_cex, point_lwd=point_lwd, ci_lwd=ci_lwd,
    ref_lwd=ref_lwd, ref_lty=ref_lty, ref_col=ref_col,
    arrow_length=arrow_length, zebra=zebra, zebra_fill=zebra_fill,
    zebra_by=zebra_by, label_width=label_width, forest_width=forest_width,
    column_gap=column_gap, panel_gap=panel_gap,
    panel_forest_ratio=panel_forest_ratio,
    label_indent=label_indent, label_wrap=label_wrap,
    width=width, height=height, dpi=dpi, file=file
  )

  if (identical(type, "multioutcome")) {
    d <- .r4vn_tf_data(data)
    specs <- .r4vn_tf_outcome_specs(outcomes, d)
    if (is.null(specs) || !length(specs)) stop("`outcomes` is required for multi-outcome mode.", call. = FALSE)
    if (is.null(predictors)) stop("`predictors` is required for multi-outcome mode.", call. = FALSE)
    panel_list <- list()
    long <- list()
    for (i in seq_along(specs)) {
      sp <- specs[[i]]
      est_i <- if (!is.null(sp$estimate)) sp$estimate else estimate[1L]
      or_i <- if (!is.null(sp$or)) isTRUE(sp$or) else isTRUE(or)
      rr_i <- if (!is.null(sp$rr)) isTRUE(sp$rr) else isTRUE(rr)
      pr_i <- if (!is.null(sp$pr)) isTRUE(sp$pr) else isTRUE(pr)
      irr_i <- if (!is.null(sp$irr)) isTRUE(sp$irr) else isTRUE(irr)
      ev_i <- if (!is.null(sp$event)) sp$event else event
      fail_i <- if (!is.null(sp$failure)) sp$failure else failure
      time_i <- if (!is.null(sp$time)) sp$time else NULL
      args <- list(
        outcome=sp$outcome, predictors=predictors, data=d,
        time=time_i, event=ev_i, failure=fail_i,
        type="regression", crude=crude, adjusted=adjusted, multi=multi,
        or=or_i, rr=rr_i, pr=pr_i, irr=irr_i, estimate=est_i,
        ci=ci, sample=sample, per=per, per_labels=per_labels,
        xmin=NULL, xmax=NULL, ticks=NULL, log=NULL, arrows=arrows,
        row_layout=row_layout, layout=layout, row_spacing=row_spacing,
        model_row_gap=model_row_gap, group_gap=group_gap, order=order,
        reference=reference, pvalue=pvalue, global_p=global_p,
        show_n=show_n, show_events=show_events,
        show_model_label=show_model_label, p_layout=p_layout,
        effect_digit=effect_digit, p_digit=p_digit,
        labels=labels, level_labels=level_labels,
        lang=lang, text=text, model_labels=model_labels,
        template=template, grid=grid, font_family=font_family,
        base_size=base_size, label_cex=label_cex, header_cex=header_cex,
        axis_cex=axis_cex, model_cex=model_cex,
        colors=colors, fills=fills, pch=pch, lty=lty,
        point_cex=point_cex, point_lwd=point_lwd, ci_lwd=ci_lwd,
        ref_lwd=ref_lwd, ref_lty=ref_lty, ref_col=ref_col,
        arrow_length=arrow_length, zebra=FALSE,
        label_width=label_width, forest_width=forest_width,
        column_gap=column_gap, label_indent=label_indent, label_wrap=label_wrap,
        show=FALSE, console=FALSE
      )
      p <- do.call(tabforest, args)
      p$panel_label <- sp$label
      panel_list[[sp$label]] <- p
      z <- p$table
      z$outcome_panel <- sp$label
      z$effect_type <- p$effect
      long[[length(long) + 1L]] <- z
    }
    keys0 <- panel_list[[1L]]$model_keys
    if (any(vapply(panel_list, function(z) !identical(z$model_keys, keys0), logical(1)))) stop("All multi-outcome panels must request the same model groups (crude/adjusted/multi).", call. = FALSE)
    rows0 <- panel_list[[1L]]$rows
    tx <- .r4vn_tf_text_full(lang, text)
    dat <- data.frame(.id=seq_len(nrow(rows0)), stringsAsFactors=FALSE)
    dat[[tx$characteristic]] <- rows0$label
    for (nm in names(panel_list)) {
      pd <- panel_list[[nm]]$data
      if (nrow(pd) != nrow(rows0)) stop("Multi-outcome panels produced incompatible row structures.", call. = FALSE)
      extra <- pd[, -1L, drop=FALSE]
      names(extra) <- paste(nm, names(extra), sep=" | ")
      dat <- data.frame(dat, extra, check.names=FALSE, stringsAsFactors=FALSE)
    }
    dat$.id <- NULL
    output <- list(data=dat, table=do.call(rbind, long), rows=rows0,
                   panels=panel_list, models=lapply(panel_list, function(z) z$models),
                   effect=vapply(panel_list, function(z) z$effect, character(1)),
                   null=vapply(panel_list, function(z) z$null, numeric(1)),
                   model_keys=keys0, settings=common_settings,
                   analysis=list(mode="multioutcome", outcomes=specs),
                   metadata=panel_list[[1L]]$metadata, call=call)
    class(output) <- c("r4vn_tabforest_multi", "r4vn_tabforest")
    if (isTRUE(console)) print(output$table, row.names=FALSE)
    if (isTRUE(show) || !is.null(file)) plot(output, file=file, width=width, height=height, dpi=dpi)
    return(invisible(output))
  }

  if (identical(type, "subgroup")) {
    d <- .r4vn_tf_data(data)
    outcome_expr <- substitute(outcome)
    outcome_name <- .r4vn_tf_name(outcome_expr, d, env, "outcome")
    time_name <- .r4vn_tf_name(substitute(time), d, env, "time", optional=TRUE)
    if (is.null(predictor)) stop("`predictor` is required for subgroup mode.", call. = FALSE)
    predictor_spec <- .r4vn_tf_spec(predictor, d, "predictor")
    if (nrow(predictor_spec) != 1L) stop("`predictor` must contain exactly one main exposure in subgroup mode.", call. = FALSE)
    subgroup_spec <- .r4vn_tf_spec(subgroup, d, "subgroup")
    effect <- .r4vn_tf_effect_type(d, outcome_name, time_name, or, rr, pr, irr, estimate)
    if (effect %in% c("OR","RR","PR") && is.null(event)) {
      lev <- .r4vn_tf_levels(d[[outcome_name]])
      if (length(lev) == 2L) event <- tail(lev,1L)
    }
    if (identical(effect,"HR") && is.null(failure)) {
      st <- d[[outcome_name]]; lev <- .r4vn_tf_levels(st)
      if (is.numeric(st) && any(st==1,na.rm=TRUE)) failure <- 1 else if (length(lev)) failure <- tail(lev,1L)
    }
    z <- .r4vn_tf_subgroup_build(d, outcome_name, time_name,
                                 predictor_spec, subgroup_spec,
                                 adjusted=adjusted, multi=multi,
                                 effect=effect, event=event, failure=failure,
                                 ci=ci, sample=sample, per=per,
                                 labels=labels, level_labels=level_labels)
    tx <- .r4vn_tf_text_full(lang, text)
    effect_key <- if (identical(effect,"Beta")) "beta" else effect
    if (is.null(effect_title)) {
      effect_title <- tx[[effect_key]]
      if (!isTRUE(all.equal(ci,.95))) effect_title <- sub("95%",paste0(format(100*ci,trim=TRUE),"%"),effect_title,fixed=TRUE)
    }
    if (is.null(label_title)) label_title <- tx$subgroup
    if (is.null(axis_title)) axis_title <- if (identical(effect,"Beta")) "Beta" else effect
    ratio <- effect %in% c("OR","RR","PR","IRR","HR")
    if (is.null(log)) log <- ratio
    common_settings$text_resolved <- tx
    common_settings$effect_title <- effect_title
    common_settings$label_title <- label_title
    common_settings$axis_title <- axis_title
    common_settings$log <- log
    publication <- .r4vn_tf_subgroup_publication(z$rows, z$estimates, tx, effect_title,
                                                  pvalue=pvalue, show_interaction_p=show_interaction_p,
                                                  effect_digit=effect_digit, p_digit=p_digit,
                                                  show_n=show_n, show_events=show_events)
    output <- list(data=publication, table=z$estimates, rows=z$rows,
                   models=z$models, interactions=z$interactions,
                   effect=effect, null=if (ratio) 1 else 0,
                   scale=if (isTRUE(log)) "log" else "linear",
                   model_keys="Subgroup", settings=common_settings,
                   analysis=list(mode="subgroup", outcome=outcome_name,
                                 time=time_name, event=event, failure=failure,
                                 predictor=predictor_spec$variable[1L],
                                 subgroups=subgroup_spec$variable,
                                 adjusted=z$adjusted_spec,
                                 common_sample=z$common_sample),
                   metadata=list(predictor=predictor_spec, subgroup=subgroup_spec), call=call)
    class(output) <- c("r4vn_tabforest_subgroup", "r4vn_tabforest")
    if (isTRUE(console)) print(output$table,row.names=FALSE)
    if (isTRUE(show) || !is.null(file)) plot(output,file=file,width=width,height=height,dpi=dpi)
    return(invisible(output))
  }

  expr <- substitute(outcome)
  obj <- if (!missing(outcome)) tryCatch(eval(expr, envir=env), error=function(e) NULL) else NULL
  object_mode <- !is.null(obj) && inherits(obj, c("r4vn_surv","r4vn_tabmulti","r4vn_stat","r4vn_result","lm","glm","coxph"))
  analysis <- list(mode=if (object_mode) "object" else "data")

  if (object_mode) {
    z <- .r4vn_tf_object(obj, select=select, ci=ci)
    estimates <- z$estimates; effect <- z$effect; models <- z$models
    model_keys <- z$model_keys; source_data <- z$data; focal_spec <- z$spec
    if (!is.null(source_data) && !is.null(focal_spec) && nrow(focal_spec) && all(focal_spec$variable %in% names(source_data))) {
      rows <- .r4vn_tf_rows(source_data, focal_spec, estimates, labels, level_labels,
                            per=NULL, per_labels=NULL, lang=lang, reference=reference)
    } else rows <- .r4vn_tf_rows_from_estimates(estimates, labels, level_labels, reference=reference)
    if (!isTRUE(reference)) estimates <- estimates[!estimates$reference,,drop=FALSE]
    analysis$source_class <- class(obj)[1L]; analysis$select <- select
  } else {
    d <- .r4vn_tf_data(data)
    outcome_name <- .r4vn_tf_name(expr,d,env,"outcome")
    time_name <- .r4vn_tf_name(substitute(time),d,env,"time",optional=TRUE)
    focal_spec <- .r4vn_tf_spec(predictors,d,"predictors")
    if (!is.null(order)) {
      order <- as.character(order)
      idx <- c(match(order,focal_spec$variable,nomatch=0L),which(!focal_spec$variable %in% order)); idx <- idx[idx>0L]
      focal_spec <- focal_spec[idx,,drop=FALSE]
    }
    effect <- .r4vn_tf_effect_type(d,outcome_name,time_name,or,rr,pr,irr,estimate)
    if (identical(effect,"HR") && is.null(time_name)) stop("HR requires `time`.",call.=FALSE)
    if (effect %in% c("OR","RR","PR") && is.null(event)) {
      lev_y <- .r4vn_tf_levels(d[[outcome_name]]); if (length(lev_y)==2L) event <- tail(lev_y,1L)
    }
    if (identical(effect,"HR") && is.null(failure)) {
      status0 <- d[[outcome_name]]; lev_s <- .r4vn_tf_levels(status0)
      if (is.numeric(status0) && any(status0==1,na.rm=TRUE)) failure <- 1 else if (length(lev_s)) failure <- tail(lev_s,1L)
    }
    if (!isTRUE(crude) && (isFALSE(adjusted)||is.null(adjusted)) && (isFALSE(multi)||is.null(multi))) stop("Request at least one model: `crude=TRUE`, `adjusted=...`, or `multi=...`.",call.=FALSE)
    z <- .r4vn_tf_build_raw(d,outcome_name,time_name,focal_spec,adjusted,multi,effect,event,failure,ci,crude,sample,per)
    estimates <- z$estimates; models <- z$models; model_keys <- z$model_keys; source_data <- d
    rows <- .r4vn_tf_rows(d,focal_spec,estimates,labels,level_labels,per,per_labels,lang,reference)
    if (!isTRUE(reference)) estimates <- estimates[!estimates$reference,,drop=FALSE]
    analysis$outcome <- outcome_name; analysis$time <- time_name; analysis$event <- event; analysis$failure <- failure
    analysis$common_sample <- z$common_sample; analysis$analysis_n <- z$analysis_n
    analysis$adjusted <- z$adjusted_spec; analysis$adjusted_all <- z$adjusted_all; analysis$multi <- z$multi_spec
  }

  effect <- unname(as.character(effect)[1L]); effect <- if (toupper(effect)=="BETA") "Beta" else toupper(effect)
  ratio <- effect %in% c("OR","RR","PR","IRR","HR"); null <- if (ratio) 1 else 0
  if (is.null(log)) log <- ratio
  if (!is.logical(log) || length(log)!=1L || is.na(log)) stop("`log` must be TRUE, FALSE, or NULL.",call.=FALSE)
  if (isTRUE(log) && any(estimates$estimate[!estimates$reference] <= 0,na.rm=TRUE)) stop("Log forest axes require positive estimates.",call.=FALSE)

  model_keys <- unique(as.character(model_keys)); display_models <- model_keys
  tx <- .r4vn_tf_text_full(lang,text)
  base_labels <- c(Crude=tx$crude,Adjusted=tx$adjusted,Multivariable=tx$multi)
  for (i in seq_along(display_models)) if (display_models[i] %in% names(base_labels)) display_models[i] <- base_labels[[display_models[i]]]
  if (!is.null(model_labels)) {
    if (is.null(names(model_labels))) display_models[seq_len(min(length(display_models),length(model_labels)))] <- as.character(model_labels)[seq_len(min(length(display_models),length(model_labels)))]
    else for (i in seq_along(model_keys)) if (model_keys[i] %in% names(model_labels)) display_models[i] <- as.character(model_labels[[model_keys[i]]])[1L]
  }
  if (is.null(label_title)) label_title <- tx$characteristic
  effect_key <- if (identical(effect,"Beta")) "beta" else effect
  if (is.null(effect_title)) {
    effect_title <- tx[[effect_key]]
    if (!isTRUE(all.equal(ci,.95))) effect_title <- sub("95%",paste0(format(100*ci,trim=TRUE),"%"),effect_title,fixed=TRUE)
  }
  if (is.null(axis_title)) axis_title <- if (identical(effect,"Beta")) "Beta" else effect

  common_settings$text_resolved <- tx
  common_settings$display_models <- stats::setNames(display_models,model_keys)
  common_settings$effect_title <- effect_title
  common_settings$label_title <- label_title
  common_settings$axis_title <- axis_title
  common_settings$log <- log

  publication_data <- .r4vn_tf_publication_data(rows,estimates,model_keys,
                                                 common_settings$display_models,effect_title,tx,
                                                 effect_digit,p_digit,pvalue,global_p,
                                                 show_n,show_events)
  output <- list(data=publication_data,table=estimates,rows=rows,models=models,
                 effect=effect,null=null,scale=if (isTRUE(log)) "log" else "linear",
                 model_keys=model_keys,settings=common_settings,analysis=analysis,
                 metadata=focal_spec,call=call)
  class(output) <- "r4vn_tabforest"
  if (isTRUE(console)) print(estimates,row.names=FALSE)
  if (isTRUE(show) || !is.null(file)) plot(output,file=file,width=width,height=height,dpi=dpi)
  invisible(output)
}


.r4vn_tf_refresh_settings <- function(x, dots = list()) {
  s <- x$settings
  for (nm in names(dots)) s[[nm]] <- dots[[nm]]
  if (is.null(s$lang)) s$lang <- "en"
  tx <- .r4vn_tf_text_full(match.arg(s$lang, c("en","vi")), s$text)
  s$text_resolved <- tx
  wording_changed <- any(c("lang","text") %in% names(dots))
  if (is.null(s$label_title) || (wording_changed && !"label_title" %in% names(dots))) {
    s$label_title <- if (inherits(x,"r4vn_tabforest_subgroup")) tx$subgroup else tx$characteristic
  }
  if (length(x$effect) == 1L) {
    key <- if (identical(x$effect,"Beta")) "beta" else x$effect
    if (is.null(s$effect_title) || (wording_changed && !"effect_title" %in% names(dots))) s$effect_title <- tx[[key]]
    if (is.null(s$axis_title)) s$axis_title <- if (identical(x$effect,"Beta")) "Beta" else x$effect
  }
  if (!is.null(x$model_keys)) {
    keys <- x$model_keys
    dm <- keys
    base <- c(Crude=tx$crude,Adjusted=tx$adjusted,Multivariable=tx$multi)
    for (i in seq_along(dm)) if (dm[i] %in% names(base)) dm[i] <- base[[dm[i]]]
    if (!is.null(s$model_labels)) {
      ml <- s$model_labels
      if (is.null(names(ml))) dm[seq_len(min(length(dm),length(ml)))] <- as.character(ml)[seq_len(min(length(dm),length(ml)))]
      else for (i in seq_along(keys)) if (keys[i] %in% names(ml)) dm[i] <- as.character(ml[[keys[i]]])[1L]
    }
    s$display_models <- stats::setNames(dm,keys)
  }
  s
}

.r4vn_tf_open_plot <- function(s, nrows, vertical_span = NULL, file = NULL,
                               width = NULL, height = NULL, dpi = NULL) {
  if (!is.null(file)) s$file <- file
  if (!is.null(width)) s$width <- width
  if (!is.null(height)) s$height <- height
  if (!is.null(dpi)) s$dpi <- dpi
  w <- suppressWarnings(as.numeric(s$width)[1L]); if (!is.finite(w) || w <= 0) w <- 12
  h <- s$height
  if (is.null(h) || !is.finite(as.numeric(h)[1L]) || as.numeric(h)[1L] <= 0) {
    base_span <- if (is.null(vertical_span) || !is.finite(vertical_span)) nrows else vertical_span
    h <- max(4.8, 2.7 + .31 * base_span)
  } else h <- as.numeric(h)[1L]
  dpi0 <- suppressWarnings(as.numeric(s$dpi)[1L]); if (!is.finite(dpi0) || dpi0 <= 0) dpi0 <- 300
  opened <- FALSE
  if (!is.null(s$file) && nzchar(as.character(s$file)[1L])) {
    .r4vn_tf_device(as.character(s$file)[1L],w,h,dpi0)
    opened <- TRUE
  }
  list(settings=s,width=w,height=h,dpi=dpi0,opened=opened)
}

.r4vn_tf_draw_axis <- function(forest_left, forest_right, axis_y, header_bottom,
                               range_obj, null, ref_col, ref_lwd, ref_lty,
                               grid = "major", axis_cex = 1, axis_title = NULL) {
  mapx <- function(v) forest_left + (range_obj$trans(v)-range_obj$tmin)/(range_obj$tmax-range_obj$tmin) * (forest_right-forest_left)
  if (grid %in% c("major","both") && length(range_obj$ticks)) {
    for (tk in range_obj$ticks) graphics::segments(mapx(tk),axis_y+.18,mapx(tk),header_bottom,col="gray90",lwd=.7)
  }
  nx <- mapx(null)
  if (is.finite(nx) && nx >= forest_left && nx <= forest_right) graphics::segments(nx,axis_y+.18,nx,header_bottom,col=ref_col,lwd=ref_lwd,lty=ref_lty)
  graphics::segments(forest_left,axis_y,forest_right,axis_y,lwd=.8)
  for (tk in range_obj$ticks) {
    xx <- mapx(tk)
    graphics::segments(xx,axis_y,xx,axis_y-.10,lwd=.8)
    lab <- if (abs(tk)>=1000 || (abs(tk)>0 && abs(tk)<.01)) format(tk,scientific=TRUE,digits=2) else formatC(tk,format="fg",digits=4,flag="#")
    lab <- sub("\\.?0+$","",lab)
    graphics::text(xx,axis_y-.25,labels=lab,cex=axis_cex,adj=c(.5,1))
  }
  if (!is.null(axis_title) && nzchar(as.character(axis_title)[1L])) graphics::text((forest_left+forest_right)/2,axis_y-.62,labels=axis_title,cex=axis_cex)
  mapx
}

.r4vn_tf_draw_ci <- function(z, py, mapx, xmin, xmax, forest_left, forest_right,
                             col, fill, pch, point_cex, point_lwd,
                             ci_lwd, lty, arrows, arrow_length) {
  if (!nrow(z) || !is.finite(z$estimate[1L]) || isTRUE(z$reference[1L])) return(FALSE)
  lo <- z$lower[1L]; hi <- z$upper[1L]; ee <- z$estimate[1L]
  if (!all(is.finite(c(lo,hi,ee)))) return(FALSE)
  lo_clip <- max(lo,xmin); hi_clip <- min(hi,xmax)
  if (hi_clip >= xmin && lo_clip <= xmax && lo_clip <= hi_clip) graphics::segments(mapx(lo_clip),py,mapx(hi_clip),py,col=col,lwd=ci_lwd,lty=lty)
  left_clip <- lo < xmin; right_clip <- hi > xmax
  if (isTRUE(arrows) && left_clip) graphics::arrows(forest_left+.012,py,forest_left,py,length=arrow_length,angle=28,code=2,col=col,lwd=ci_lwd)
  if (isTRUE(arrows) && right_clip) graphics::arrows(forest_right-.012,py,forest_right,py,length=arrow_length,angle=28,code=2,col=col,lwd=ci_lwd)
  if (ee >= xmin && ee <= xmax) graphics::points(mapx(ee),py,pch=pch,cex=point_cex,col=col,bg=fill,lwd=point_lwd)
  else if (isTRUE(arrows)) {
    if (ee < xmin) graphics::arrows(forest_left+.025,py,forest_left,py,length=arrow_length,angle=28,code=2,col=col,lwd=ci_lwd*1.15)
    if (ee > xmax) graphics::arrows(forest_right-.025,py,forest_right,py,length=arrow_length,angle=28,code=2,col=col,lwd=ci_lwd*1.15)
  }
  left_clip || right_clip || ee < xmin || ee > xmax
}

#' Plot an R4VN regression forest object
#'
#' Re-draws a `tabforest()` object without refitting models. All graphical
#' settings may be overridden here, which is useful for trying different journal
#' layouts after the analysis is finalized.
#'
#' @param x An `r4vn_tabforest` object.
#' @param ... Any plotting option accepted by `tabforest()`.
#' @param file,width,height,dpi Optional output overrides.
#' @return `x`, invisibly.
#' @export
plot.r4vn_tabforest <- function(x, ..., file=NULL, width=NULL, height=NULL, dpi=NULL) {
  if (inherits(x,"r4vn_tabforest_multi")) return(plot.r4vn_tabforest_multi(x,...,file=file,width=width,height=height,dpi=dpi))
  if (inherits(x,"r4vn_tabforest_subgroup")) return(plot.r4vn_tabforest_subgroup(x,...,file=file,width=width,height=height,dpi=dpi))
  if (!inherits(x,"r4vn_tabforest")) stop("`x` must be an r4vn_tabforest object.",call.=FALSE)
  s <- .r4vn_tf_refresh_settings(x,list(...))
  if (!is.null(file)) s$file <- file; if (!is.null(width)) s$width <- width; if (!is.null(height)) s$height <- height; if (!is.null(dpi)) s$dpi <- dpi
  tx <- s$text_resolved; est <- x$table; rows <- x$rows; keys <- x$model_keys
  nmodels <- length(keys); if (!nmodels) stop("The forest contains no model columns.",call.=FALSE)
  ratio <- x$effect %in% c("OR","RR","PR","IRR","HR")
  use_log <- isTRUE(s$log)
  rng <- .r4vn_tf_axis_range(est,x$null,s$xmin,s$xmax,s$ticks,use_log,ratio)

  dm <- s$display_models; if (is.null(dm)) dm <- stats::setNames(keys,keys)
  row_layout_plot <- s$row_layout
  if (is.null(row_layout_plot) || identical(row_layout_plot,"auto")) row_layout_plot <- if (nmodels > 1L) "modelrows" else "compact"
  rows_plot <- .r4vn_tf_expand_rows(rows,keys,dm,row_layout_plot,s$show_model_label)
  row_y <- .r4vn_tf_row_positions(rows_plot,s$row_spacing,s$model_row_gap,s$group_gap)
  nrows <- nrow(rows_plot); y_top <- max(row_y)+1.7; y_bottom <- -.15
  vertical_span <- max(row_y)-min(row_y)+5
  op <- .r4vn_tf_open_plot(s,nrows,vertical_span,file,width,height,dpi); s <- op$settings; opened <- op$opened
  oldpar <- graphics::par(no.readonly=TRUE)
  on.exit({try(graphics::par(oldpar),silent=TRUE); if (opened) try(grDevices::dev.off(),silent=TRUE)},add=TRUE)
  family <- if (is.null(s$font_family)) "" else as.character(s$font_family)[1L]
  graphics::par(mar=c(1.6,1,2.8,1),xaxs="i",yaxs="i",family=family,ps=as.numeric(s$base_size))
  graphics::plot.new(); graphics::plot.window(xlim=c(0,1),ylim=c(y_bottom,y_top),xaxs="i",yaxs="i")
  graphics::rect(0,y_bottom,1,y_top,col="white",border=NA)
  if (isTRUE(s$zebra)) .r4vn_tf_draw_zebra(rows_plot,row_y,s$zebra_fill,s$zebra_by)

  lw <- as.numeric(s$label_width); fw <- as.numeric(s$forest_width)
  if (!is.finite(lw) || lw<=.12 || lw>=.60) stop("`label_width` should be between 0.12 and 0.60.",call.=FALSE)
  if (!is.finite(fw) || fw<=.12 || fw>=.60) stop("`forest_width` should be between 0.12 and 0.60.",call.=FALSE)
  forest_left <- lw; forest_right <- forest_left+fw
  if (forest_right>=.86) stop("`label_width + forest_width` leaves too little room for numbers.",call.=FALSE)
  num_left <- forest_right + as.numeric(s$column_gap); num_right <- .995
  header_y <- max(row_y)+1.0; header_bottom <- header_y-.60; axis_y <- .58
  graphics::text(.008,header_y,labels=s$label_title,adj=c(0,.5),font=2,cex=s$header_cex)
  graphics::text((forest_left+forest_right)/2,header_y,labels=s$effect_title,font=2,cex=s$header_cex)

  modelrows <- identical(row_layout_plot,"modelrows")
  if (modelrows) {
    if (identical(s$p_layout,"column") && isTRUE(s$pvalue)) {
      eff_right <- num_left + (num_right-num_left)*.76
      p_left <- eff_right
      graphics::text((num_left+eff_right)/2,header_y,labels=s$effect_title,font=2,cex=s$header_cex*.88)
      graphics::text((p_left+num_right)/2,header_y,labels=tx$p,font=2,cex=s$header_cex*.88)
    } else {
      eff_right <- num_right; p_left <- num_right
      graphics::text((num_left+num_right)/2,header_y,labels=s$effect_title,font=2,cex=s$header_cex*.88)
    }
  } else {
    model_width <- (num_right-num_left)/nmodels
    for (j in seq_along(keys)) {
      a <- num_left+(j-1)*model_width; b <- a+model_width
      graphics::text((a+b)/2,header_y+.20,labels=as.character(dm[[keys[j]]]),font=2,cex=s$header_cex*.88)
      graphics::text((a+b)/2,header_y-.22,labels=s$effect_title,font=2,cex=s$header_cex*.76)
    }
  }
  graphics::segments(.005,header_bottom,.995,header_bottom,lwd=1)
  mapx <- .r4vn_tf_draw_axis(forest_left,forest_right,axis_y,header_bottom,rng,x$null,
                             s$ref_col,s$ref_lwd,s$ref_lty,s$grid,s$axis_cex,s$axis_title)

  cols <- .r4vn_tf_resolve_named(s$colors,keys,"black")
  fills <- .r4vn_tf_resolve_named(s$fills,keys,NA)
  pchs <- as.numeric(.r4vn_tf_resolve_named(s$pch,keys,c(1,19,15,17)[seq_len(nmodels)]))
  ltys <- as.numeric(.r4vn_tf_resolve_named(s$lty,keys,1))
  pcex <- as.numeric(.r4vn_tf_resolve_named(s$point_cex,keys,1.15))
  plwd <- as.numeric(.r4vn_tf_resolve_named(s$point_lwd,keys,1))
  clwd <- as.numeric(.r4vn_tf_resolve_named(s$ci_lwd,keys,1.2))
  dodge <- if (nmodels<=1L) 0 else if (identical(s$layout,"stack")) .20*as.numeric(s$row_spacing) else .11*as.numeric(s$row_spacing)
  offsets <- if (nmodels<=1L) 0 else seq(dodge,-dodge,length.out=nmodels)
  clipped_any <- FALSE

  for (i in seq_len(nrows)) {
    r <- rows_plot[i,,drop=FALSE]; yy <- row_y[i]
    if (identical(r$row_type,"header")) {
      graphics::text(.008,yy,labels=.r4vn_tf_wrap(r$label,s$label_wrap),adj=c(0,.5),font=2,cex=s$label_cex)
      if (isTRUE(s$global_p)) {
        pieces <- character()
        for (k in keys) {
          z <- est[est$variable==r$variable & est$model==k,,drop=FALSE]
          gp <- if (nrow(z)) z$global_p[is.finite(z$global_p)] else numeric()
          if (length(gp)) pieces <- c(pieces,paste0(as.character(dm[[k]])," ",tx$p,"=",.r4vn_tf_fmt_p(gp[1L],s$p_digit)))
        }
        if (length(pieces)) graphics::text(num_left,yy,labels=paste(pieces,collapse="; "),adj=c(0,.5),cex=s$label_cex*.82)
      }
      next
    }

    lx <- if (identical(r$row_type,"level")) .008+as.numeric(s$label_indent) else .008
    label_here <- r$label
    if (modelrows && i>1L && isTRUE(rows_plot$.item_id[i]==rows_plot$.item_id[i-1L])) label_here <- ""
    graphics::text(lx,yy,labels=.r4vn_tf_wrap(label_here,s$label_wrap),adj=c(0,.5),cex=s$label_cex)
    if (modelrows && !is.na(r$.display_model) && nzchar(r$.display_model_label)) {
      graphics::text(forest_left-.012,yy,labels=r$.display_model_label,adj=c(1,.5),cex=s$model_cex,font=3)
    }

    draw_keys <- if (modelrows) r$.display_model else keys
    for (k in draw_keys) {
      if (is.na(k) || !nzchar(k)) next
      j <- match(k,keys)
      if (identical(r$row_type,"level")) z <- est[est$variable==r$variable & est$level==r$level & est$model==k,,drop=FALSE]
      else z <- est[est$variable==r$variable & (is.na(est$level)|est$level=="") & est$model==k,,drop=FALSE]
      if (!nrow(z)) next
      ref <- isTRUE(z$reference[1L]); py <- if (modelrows) yy else yy+offsets[j]
      nt <- .r4vn_tf_numeric_text(z,ref,tx,s$effect_digit,s$p_digit,s$pvalue,s$show_n,s$show_events,s$p_layout)
      if (modelrows) {
        if (identical(s$p_layout,"column") && isTRUE(s$pvalue)) {
          graphics::text(num_left,yy,labels=nt$effect,adj=c(0,.5),cex=s$label_cex*.84)
          if (nzchar(nt$p)) graphics::text((p_left+num_right)/2,yy,labels=nt$p,cex=s$label_cex*.84)
        } else graphics::text(num_left,yy,labels=nt$effect,adj=c(0,.5),cex=s$label_cex*.84)
      } else {
        model_width <- (num_right-num_left)/nmodels; a <- num_left+(j-1)*model_width; b <- a+model_width
        if (identical(s$p_layout,"column") && isTRUE(s$pvalue)) {
          split <- a+(b-a)*.77
          graphics::text(a+.003,py,labels=nt$effect,adj=c(0,.5),cex=s$label_cex*.78)
          if (nzchar(nt$p)) graphics::text((split+b)/2,py,labels=nt$p,cex=s$label_cex*.78)
        } else graphics::text((a+b)/2,py,labels=nt$effect,cex=s$label_cex*.78)
      }
      clipped_any <- .r4vn_tf_draw_ci(z,py,mapx,rng$xmin,rng$xmax,forest_left,forest_right,
                                       cols[[k]],fills[[k]],pchs[j],pcex[j],plwd[j],clwd[j],ltys[j],
                                       s$arrows,s$arrow_length) || clipped_any
    }
  }

  if (!is.null(s$title) && nzchar(as.character(s$title)[1L])) graphics::title(main=as.character(s$title)[1L],cex.main=1.10)
  if (!is.null(s$subtitle) && nzchar(as.character(s$subtitle)[1L])) graphics::mtext(as.character(s$subtitle)[1L],side=3,line=.15,cex=.90)
  notes <- character()
  if (!is.null(s$caption) && nzchar(as.character(s$caption)[1L])) notes <- c(notes,as.character(s$caption)[1L])
  if (is.character(s$note) && length(s$note) && nzchar(s$note[1L])) notes <- c(notes,s$note[1L]) else if (isTRUE(s$note) && clipped_any) notes <- c(notes,tx$arrow_note)
  if (length(notes)) graphics::mtext(paste(notes,collapse="\n"),side=1,line=.2,adj=0,cex=.78)
  invisible(x)
}

#' @export
plot.r4vn_tabforest_multi <- function(x, ..., file=NULL, width=NULL, height=NULL, dpi=NULL) {
  if (!inherits(x,"r4vn_tabforest_multi")) stop("`x` must be an r4vn_tabforest_multi object.",call.=FALSE)
  s <- .r4vn_tf_refresh_settings(x,list(...))
  if (!is.null(file)) s$file<-file; if (!is.null(width)) s$width<-width; if (!is.null(height)) s$height<-height; if (!is.null(dpi)) s$dpi<-dpi
  tx <- s$text_resolved; panels <- x$panels; pnames <- names(panels); np <- length(panels); keys <- x$model_keys; nmodels <- length(keys)
  dm <- s$display_models; if (is.null(dm)) dm <- stats::setNames(keys,keys)
  row_layout_plot <- s$row_layout
  if (is.null(row_layout_plot) || identical(row_layout_plot,"auto")) row_layout_plot <- if (nmodels > 1L) "modelrows" else "compact"
  rows_plot <- .r4vn_tf_expand_rows(x$rows,keys,dm,row_layout_plot,s$show_model_label)
  row_y <- .r4vn_tf_row_positions(rows_plot,s$row_spacing,s$model_row_gap,s$group_gap)
  nrows <- nrow(rows_plot); y_top <- max(row_y)+1.85; y_bottom <- -.18
  op <- .r4vn_tf_open_plot(s,nrows,max(row_y)-min(row_y)+5,file,width,height,dpi); s<-op$settings; opened<-op$opened
  oldpar <- graphics::par(no.readonly=TRUE)
  on.exit({try(graphics::par(oldpar),silent=TRUE);if(opened)try(grDevices::dev.off(),silent=TRUE)},add=TRUE)
  family <- if(is.null(s$font_family))"" else as.character(s$font_family)[1L]
  graphics::par(mar=c(1.8,1,3,1),xaxs="i",yaxs="i",family=family,ps=as.numeric(s$base_size))
  graphics::plot.new();graphics::plot.window(xlim=c(0,1),ylim=c(y_bottom,y_top),xaxs="i",yaxs="i")
  graphics::rect(0,y_bottom,1,y_top,col="white",border=NA)
  if(isTRUE(s$zebra)) .r4vn_tf_draw_zebra(rows_plot,row_y,s$zebra_fill,s$zebra_by)

  lw <- as.numeric(s$label_width); if(!is.finite(lw)||lw<=.12||lw>=.55) stop("For multi-outcome plots, `label_width` should be between 0.12 and 0.55.",call.=FALSE)
  pg <- as.numeric(s$panel_gap); if(!is.finite(pg)||pg<0)pg<-.012
  usable <- 1-lw-pg*(np-1L)-.01; if(usable<=.25) stop("Too many panels for the current `width`/`label_width`.",call.=FALSE)
  pw <- usable/np; fr <- as.numeric(s$panel_forest_ratio); if(!is.finite(fr)||fr<.35||fr>.80)fr<-.58
  header_y <- max(row_y)+1.08; header_bottom <- header_y-.63; axis_y <- .60
  graphics::text(.008,header_y,labels=if(is.null(s$label_title))tx$characteristic else s$label_title,adj=c(0,.5),font=2,cex=s$header_cex)
  graphics::segments(.005,header_bottom,.995,header_bottom,lwd=1)

  cols <- .r4vn_tf_resolve_named(s$colors,keys,"black"); fills <- .r4vn_tf_resolve_named(s$fills,keys,NA)
  pchs <- as.numeric(.r4vn_tf_resolve_named(s$pch,keys,c(1,19,15,17)[seq_len(nmodels)])); ltys <- as.numeric(.r4vn_tf_resolve_named(s$lty,keys,1))
  pcex <- as.numeric(.r4vn_tf_resolve_named(s$point_cex,keys,1.15)); plwd <- as.numeric(.r4vn_tf_resolve_named(s$point_lwd,keys,1)); clwd <- as.numeric(.r4vn_tf_resolve_named(s$ci_lwd,keys,1.2))
  dodge <- if(nmodels<=1L)0 else if(identical(s$layout,"stack")) .20*as.numeric(s$row_spacing) else .11*as.numeric(s$row_spacing)
  offsets <- if(nmodels<=1L)0 else seq(dodge,-dodge,length.out=nmodels)
  panel_geom <- vector("list",np); clipped_any <- FALSE

  for(q in seq_along(panels)) {
    p <- panels[[q]]; nm <- pnames[q]
    left <- lw+(q-1L)*(pw+pg); right <- left+pw
    fleft <- left+.015*pw; fright <- left+fr*pw
    nleft <- fright+.025*pw; nright <- right-.010*pw
    effect <- p$effect; ratio <- effect %in% c("OR","RR","PR","IRR","HR")
    logq <- .r4vn_tf_panel_value(s$log,nm,q,ratio); if(is.null(logq))logq<-ratio
    xminq <- .r4vn_tf_panel_value(s$xmin,nm,q,NULL); xmaxq <- .r4vn_tf_panel_value(s$xmax,nm,q,NULL)
    ticksq <- .r4vn_tf_panel_value(s$ticks,nm,q,NULL)
    rng <- .r4vn_tf_axis_range(p$table,p$null,xminq,xmaxq,ticksq,isTRUE(logq),ratio)
    graphics::text((left+right)/2,header_y+.28,labels=nm,font=2,cex=s$header_cex*.90)
    effect_key <- if (identical(effect,"Beta")) "beta" else effect
    etitle <- if (!is.null(s$effect_title)) s$effect_title else tx[[effect_key]]
    graphics::text((fleft+fright)/2,header_y-.18,labels=etitle,font=2,cex=s$header_cex*.72)
    axis_q <- if (!is.null(s$axis_title)) s$axis_title else if (identical(effect,"Beta")) "Beta" else effect
    mapx <- .r4vn_tf_draw_axis(fleft,fright,axis_y,header_bottom,rng,p$null,s$ref_col,s$ref_lwd,s$ref_lty,s$grid,s$axis_cex,axis_q)
    if(q>1L) graphics::segments(left-pg/2,axis_y-.35,left-pg/2,header_y+.55,col="gray85",lwd=.7)
    panel_geom[[q]] <- list(fleft=fleft,fright=fright,nleft=nleft,nright=nright,rng=rng,mapx=mapx,panel=p)
  }

  modelrows <- identical(row_layout_plot,"modelrows")
  for(i in seq_len(nrows)) {
    r <- rows_plot[i,,drop=FALSE]; yy<-row_y[i]
    if(identical(r$row_type,"header")) {
      graphics::text(.008,yy,labels=.r4vn_tf_wrap(r$label,s$label_wrap),adj=c(0,.5),font=2,cex=s$label_cex)
      next
    }
    lx <- if(identical(r$row_type,"level")) .008+as.numeric(s$label_indent) else .008
    lab <- r$label; if(modelrows && i>1L && isTRUE(rows_plot$.item_id[i]==rows_plot$.item_id[i-1L])) lab<-""
    graphics::text(lx,yy,labels=.r4vn_tf_wrap(lab,s$label_wrap),adj=c(0,.5),cex=s$label_cex)
    if(modelrows && !is.na(r$.display_model) && nzchar(r$.display_model_label)) graphics::text(lw-.010,yy,labels=r$.display_model_label,adj=c(1,.5),cex=s$model_cex,font=3)
    draw_keys <- if(modelrows) r$.display_model else keys
    for(q in seq_along(panels)) {
      pgm <- panel_geom[[q]]; p <- pgm$panel; est <- p$table
      for(k in draw_keys) {
        if(is.na(k)||!nzchar(k))next; j<-match(k,keys)
        if(identical(r$row_type,"level")) z<-est[est$variable==r$variable & est$level==r$level & est$model==k,,drop=FALSE]
        else z<-est[est$variable==r$variable & (is.na(est$level)|est$level=="") & est$model==k,,drop=FALSE]
        if(!nrow(z))next
        py <- if(modelrows)yy else yy+offsets[j]; ref<-isTRUE(z$reference[1L])
        nt <- .r4vn_tf_numeric_text(z,ref,tx,s$effect_digit,s$p_digit,s$pvalue,s$show_n,s$show_events,"inline")
        if(!modelrows && nmodels>1L) nt$effect <- paste0(as.character(dm[[k]]),": ",nt$effect)
        graphics::text(pgm$nleft,py,labels=nt$effect,adj=c(0,.5),cex=s$label_cex*.70)
        clipped_any <- .r4vn_tf_draw_ci(z,py,pgm$mapx,pgm$rng$xmin,pgm$rng$xmax,pgm$fleft,pgm$fright,
                                         cols[[k]],fills[[k]],pchs[j],pcex[j],plwd[j],clwd[j],ltys[j],s$arrows,s$arrow_length) || clipped_any
      }
    }
  }
  if(!is.null(s$title)&&nzchar(as.character(s$title)[1L]))graphics::title(main=as.character(s$title)[1L],cex.main=1.08)
  if(!is.null(s$subtitle)&&nzchar(as.character(s$subtitle)[1L]))graphics::mtext(as.character(s$subtitle)[1L],side=3,line=.15,cex=.90)
  notes<-character();if(!is.null(s$caption)&&nzchar(as.character(s$caption)[1L]))notes<-c(notes,as.character(s$caption)[1L])
  if(is.character(s$note)&&length(s$note)&&nzchar(s$note[1L]))notes<-c(notes,s$note[1L]) else if(isTRUE(s$note)&&clipped_any)notes<-c(notes,tx$arrow_note)
  if(length(notes))graphics::mtext(paste(notes,collapse="\n"),side=1,line=.25,adj=0,cex=.75)
  invisible(x)
}

#' @export
plot.r4vn_tabforest_subgroup <- function(x, ..., file=NULL, width=NULL, height=NULL, dpi=NULL) {
  if(!inherits(x,"r4vn_tabforest_subgroup"))stop("`x` must be an r4vn_tabforest_subgroup object.",call.=FALSE)
  s <- .r4vn_tf_refresh_settings(x,list(...)); if(!is.null(file))s$file<-file;if(!is.null(width))s$width<-width;if(!is.null(height))s$height<-height;if(!is.null(dpi))s$dpi<-dpi
  tx<-s$text_resolved;est<-x$table;rows<-x$rows;ratio<-x$effect %in% c("OR","RR","PR","IRR","HR");use_log<-isTRUE(s$log)
  rng<-.r4vn_tf_axis_range(est,x$null,s$xmin,s$xmax,s$ticks,use_log,ratio)
  rows_plot<-rows;rows_plot$.item_id<-seq_len(nrow(rows));rows_plot$.row_group<-if(nrow(rows)) cumsum(c(TRUE,rows$variable[-1L]!=rows$variable[-nrow(rows)])) else integer()
  row_y<-.r4vn_tf_row_positions(rows_plot,s$row_spacing,s$model_row_gap,s$group_gap);nrows<-nrow(rows_plot);y_top<-max(row_y)+1.7;y_bottom<--.15
  op<-.r4vn_tf_open_plot(s,nrows,max(row_y)-min(row_y)+5,file,width,height,dpi);s<-op$settings;opened<-op$opened
  oldpar<-graphics::par(no.readonly=TRUE);on.exit({try(graphics::par(oldpar),silent=TRUE);if(opened)try(grDevices::dev.off(),silent=TRUE)},add=TRUE)
  family<-if(is.null(s$font_family))"" else as.character(s$font_family)[1L]
  graphics::par(mar=c(1.6,1,2.8,1),xaxs="i",yaxs="i",family=family,ps=as.numeric(s$base_size));graphics::plot.new();graphics::plot.window(xlim=c(0,1),ylim=c(y_bottom,y_top),xaxs="i",yaxs="i")
  graphics::rect(0,y_bottom,1,y_top,col="white",border=NA);if(isTRUE(s$zebra)).r4vn_tf_draw_zebra(rows_plot,row_y,s$zebra_fill,s$zebra_by)
  lw<-as.numeric(s$label_width);fw<-as.numeric(s$forest_width);forest_left<-lw;forest_right<-forest_left+fw
  if(forest_right>=.72)stop("For subgroup plots, reduce `label_width` or `forest_width` to leave room for p-value columns.",call.=FALSE)
  rem_left<-forest_right+as.numeric(s$column_gap);rem_right<-.995
  showp<-isTRUE(s$pvalue);showip<-isTRUE(s$show_interaction_p)
  if(showp&&showip){eff_right<-rem_left+(rem_right-rem_left)*.58;p_right<-rem_left+(rem_right-rem_left)*.78;ip_right<-rem_right}
  else if(showp){eff_right<-rem_left+(rem_right-rem_left)*.72;p_right<-rem_right;ip_right<-rem_right}
  else if(showip){eff_right<-rem_left+(rem_right-rem_left)*.70;p_right<-eff_right;ip_right<-rem_right}
  else {eff_right<-rem_right;p_right<-rem_right;ip_right<-rem_right}
  header_y<-max(row_y)+1.0;header_bottom<-header_y-.60;axis_y<-.58
  graphics::text(.008,header_y,labels=s$label_title,adj=c(0,.5),font=2,cex=s$header_cex)
  graphics::text((forest_left+forest_right)/2,header_y,labels=s$effect_title,font=2,cex=s$header_cex)
  graphics::text((rem_left+eff_right)/2,header_y,labels=s$effect_title,font=2,cex=s$header_cex*.82)
  if(showp)graphics::text((eff_right+p_right)/2,header_y,labels=tx$p,font=2,cex=s$header_cex*.82)
  if(showip)graphics::text((max(p_right,eff_right)+ip_right)/2,header_y,labels=tx$interaction_p,font=2,cex=s$header_cex*.78)
  graphics::segments(.005,header_bottom,.995,header_bottom,lwd=1)
  mapx<-.r4vn_tf_draw_axis(forest_left,forest_right,axis_y,header_bottom,rng,x$null,s$ref_col,s$ref_lwd,s$ref_lty,s$grid,s$axis_cex,s$axis_title)
  col<-.r4vn_tf_resolve_named(s$colors,"Subgroup","black")[[1L]];fill<-.r4vn_tf_resolve_named(s$fills,"Subgroup",NA)[[1L]];pch0<-as.numeric(.r4vn_tf_resolve_named(s$pch,"Subgroup",15)[[1L]])
  lty0<-as.numeric(.r4vn_tf_resolve_named(s$lty,"Subgroup",1)[[1L]]);pcex<-as.numeric(.r4vn_tf_resolve_named(s$point_cex,"Subgroup",1.1)[[1L]]);plwd<-as.numeric(.r4vn_tf_resolve_named(s$point_lwd,"Subgroup",1)[[1L]]);clwd<-as.numeric(.r4vn_tf_resolve_named(s$ci_lwd,"Subgroup",1.2)[[1L]])
  clipped_any<-FALSE
  for(i in seq_len(nrows)){
    r<-rows_plot[i,,drop=FALSE];yy<-row_y[i]
    if(identical(r$row_type,"header")){
      graphics::text(.008,yy,labels=.r4vn_tf_wrap(r$label,s$label_wrap),adj=c(0,.5),font=2,cex=s$label_cex)
      if(showip){z<-est[est$subgroup_variable==r$variable,,drop=FALSE];ip<-z$interaction_p[is.finite(z$interaction_p)];if(length(ip))graphics::text((max(p_right,eff_right)+ip_right)/2,yy,labels=.r4vn_tf_fmt_p(ip[1L],s$p_digit),cex=s$label_cex*.86)}
      next
    }
    graphics::text(.008+as.numeric(s$label_indent),yy,labels=.r4vn_tf_wrap(r$label,s$label_wrap),adj=c(0,.5),cex=s$label_cex)
    z<-est[est$subgroup_variable==r$variable & est$subgroup_level==r$level,,drop=FALSE]
    if(!nrow(z))next
    nt<-.r4vn_tf_numeric_text(z,FALSE,tx,s$effect_digit,s$p_digit,FALSE,s$show_n,s$show_events,"column")
    graphics::text(rem_left,yy,labels=nt$effect,adj=c(0,.5),cex=s$label_cex*.82)
    if(showp&&is.finite(z$p[1L]))graphics::text((eff_right+p_right)/2,yy,labels=.r4vn_tf_fmt_p(z$p[1L],s$p_digit),cex=s$label_cex*.84)
    clipped_any<-.r4vn_tf_draw_ci(z,yy,mapx,rng$xmin,rng$xmax,forest_left,forest_right,col,fill,pch0,pcex,plwd,clwd,lty0,s$arrows,s$arrow_length)||clipped_any
  }
  if(!is.null(s$title)&&nzchar(as.character(s$title)[1L]))graphics::title(main=as.character(s$title)[1L],cex.main=1.08)
  if(!is.null(s$subtitle)&&nzchar(as.character(s$subtitle)[1L]))graphics::mtext(as.character(s$subtitle)[1L],side=3,line=.15,cex=.90)
  notes<-character();if(!is.null(s$caption)&&nzchar(as.character(s$caption)[1L]))notes<-c(notes,as.character(s$caption)[1L]);if(is.character(s$note)&&length(s$note)&&nzchar(s$note[1L]))notes<-c(notes,s$note[1L]) else if(isTRUE(s$note)&&clipped_any)notes<-c(notes,tx$arrow_note);if(length(notes))graphics::mtext(paste(notes,collapse="\n"),side=1,line=.2,adj=0,cex=.78)
  invisible(x)
}

#' @export
print.r4vn_tabforest <- function(x, ...) { plot(x,...); invisible(x) }

#' @export
as.data.frame.r4vn_tabforest <- function(x, ...) x$data