R/models.R

In R4VN: Health Data Analysis and Publication-Ready Reporting

# R4VN correlation and regression models
#
# corr(), regress(), logistic(), poisson(), and lrtest() are public functions.
# Shared covariance and display helpers remain in statistics-utils.R.
#
# Compact regression syntax is implemented locally in this file so table-oriented
# vars() and the general statistical utility layer remain unchanged.

# ============================================================================
# Compact model syntax helpers
# ============================================================================

.r4vn_mx_deparse1 <- function(x) {
  paste(deparse(x, width.cutoff = 500L), collapse = "")
}

.r4vn_mx_directive_store <- function() {
  new.env(parent = emptyenv())
}


# Return the base variable name used by compact categorical prefixes.
.r4vn_mx_prefixed_base <- function(x) {
  if (!is.symbol(x)) return(NULL)
  nm <- as.character(x)

  if (grepl("^ib[0-9]+\\.", nm)) {
    return(sub("^ib[0-9]+\\.", "", nm))
  }
  if (grepl("^b[1-9][0-9]*\\.", nm)) {
    return(sub("^b[1-9][0-9]*\\.", "", nm))
  }
  if (startsWith(nm, "i.") && nchar(nm) > 2L) {
    return(sub("^i\\.", "", nm))
  }

  NULL
}

# Collect categorical variables appearing in a compact expression.
.r4vn_mx_collect_prefixed <- function(x) {
  if (is.symbol(x)) {
    z <- .r4vn_mx_prefixed_base(x)
    return(if (is.null(z)) character() else z)
  }

  if (!is.call(x) || length(x) <= 1L) return(character())

  unique(unlist(
    lapply(as.list(x)[-1L], .r4vn_mx_collect_prefixed),
    use.names = FALSE
  ))
}

# For starred categorical interactions, R's `*` already includes the main
# effects. This helper identifies those categorical main effects so redundant
# plain variables can be removed before the formula is built.
.r4vn_mx_categorical_main_bases <- function(x) {
  if (is.symbol(x)) {
    z <- .r4vn_mx_prefixed_base(x)
    return(if (is.null(z)) character() else z)
  }

  if (!is.call(x) || length(x) <= 1L) return(character())

  hd <- x[[1L]]

  if (is.symbol(hd) && identical(as.character(hd), "*")) {
    return(.r4vn_mx_collect_prefixed(x))
  }

  if (is.symbol(hd) && identical(as.character(hd), "+")) {
    return(unique(unlist(
      lapply(as.list(x)[-1L], .r4vn_mx_categorical_main_bases),
      use.names = FALSE
    )))
  }

  character()
}

.r4vn_mx_drop_redundant_plain <- function(exprs) {
  if (!length(exprs)) return(exprs)

  categorical_main <- unique(unlist(
    lapply(exprs, .r4vn_mx_categorical_main_bases),
    use.names = FALSE
  ))

  if (!length(categorical_main)) return(exprs)

  keep <- vapply(
    exprs,
    function(x) {
      if (!is.symbol(x)) return(TRUE)

      nm <- as.character(x)

      # Explicit compact declarations are not redundant plain terms.
      if (grepl("^(ib[0-9]+|b[1-9][0-9]*|i|c)\\.", nm)) {
        return(TRUE)
      }

      !nm %in% categorical_main
    },
    logical(1)
  )

  exprs[keep]
}

# Backward-compatibility alias for older model-syntax tests and integrations.
#
# R4VN now applies reference categories directly to the model data so the
# fitted model frame retains ordinary variable names such as `occupation`.
# Older code looked for an implementation-detail column whose name contained
# `r4vn_factor_index`. We add a harmless alias to fit$model after fitting.
# It does not enter the formula, design matrix, likelihood, or coefficients.
.r4vn_mx_add_reference_aliases <- function(fit, formula) {
  if (is.null(fit$model) || !is.data.frame(fit$model)) return(fit)

  directives <- attr(formula, "r4vn_model_directives")
  if (is.null(directives) || !length(directives)) return(fit)

  for (nm in names(directives)) {
    d <- directives[[nm]]

    if (is.null(d$ref) || !nm %in% names(fit$model)) next

    ref_txt <- as.character(d$ref)[1L]
    hybrid <- identical(d$ref_mode, "hybrid")

    alias <- paste0(
      ".r4vn_factor_index(",
      nm, ", ", ref_txt, ", ",
      if (hybrid) "TRUE" else "FALSE",
      ")"
    )

    if (!alias %in% names(fit$model)) {
      fit$model[[alias]] <- fit$model[[nm]]
    }
  }

  fit
}

.r4vn_mx_add_directive <- function(store, variable,
                                   type = c("factor", "continuous"),
                                   ref = NULL,
                                   ref_mode = NULL) {
  type <- match.arg(type)
  if (!nzchar(variable)) {
    stop("A model prefix must be followed by a variable name.", call. = FALSE)
  }

  incoming <- list(
    variable = variable,
    type = type,
    ref = ref,
    ref_mode = ref_mode
  )

  if (!exists(variable, envir = store, inherits = FALSE)) {
    assign(variable, incoming, envir = store)
    return(invisible(NULL))
  }

  current <- get(variable, envir = store, inherits = FALSE)

  if (!identical(current$type, incoming$type)) {
    stop(
      sprintf(
        "Variable `%s` has conflicting model prefixes. Do not use it as both continuous and categorical.",
        variable
      ),
      call. = FALSE
    )
  }

  if (identical(type, "continuous")) {
    return(invisible(NULL))
  }

  if (!is.null(current$ref) && is.null(incoming$ref)) {
    return(invisible(NULL))
  }

  if (is.null(current$ref) && !is.null(incoming$ref)) {
    assign(variable, incoming, envir = store)
    return(invisible(NULL))
  }

  if (!is.null(current$ref) && !is.null(incoming$ref)) {
    same_ref <- identical(as.character(current$ref), as.character(incoming$ref))
    same_mode <- identical(current$ref_mode, incoming$ref_mode)
    if (!same_ref || !same_mode) {
      stop(
        sprintf(
          "Variable `%s` has conflicting reference-category declarations in the same model.",
          variable
        ),
        call. = FALSE
      )
    }
  }

  invisible(NULL)
}

.r4vn_mx_collect_directives <- function(store) {
  nm <- ls(envir = store, all.names = TRUE)
  if (!length(nm)) return(list())
  out <- lapply(nm, function(z) get(z, envir = store, inherits = FALSE))
  names(out) <- nm
  out
}

.r4vn_mx_translate_expr <- function(x, store) {
  if (is.symbol(x)) {
    txt <- as.character(x)

    # ib#.var: prefer the literal value/level; if absent, a positive integer
    # can fall back to the corresponding factor-level position.
    if (grepl("^ib[0-9]+\\.", txt)) {
      ref <- sub("^ib([0-9]+)\\..*$", "\\1", txt)
      variable <- sub("^ib[0-9]+\\.", "", txt)
      .r4vn_mx_add_directive(
        store, variable, "factor",
        ref = ref, ref_mode = "hybrid"
      )
      return(as.name(variable))
    }

    # Existing R4VN bN.var convention: Nth factor level is the reference.
    if (grepl("^b[1-9][0-9]*\\.", txt)) {
      ref <- as.integer(sub("^b([1-9][0-9]*)\\..*$", "\\1", txt))
      variable <- sub("^b[1-9][0-9]*\\.", "", txt)
      .r4vn_mx_add_directive(
        store, variable, "factor",
        ref = ref, ref_mode = "index"
      )
      return(as.name(variable))
    }

    # i.var: categorical.
    if (startsWith(txt, "i.") && nchar(txt) > 2L) {
      variable <- substring(txt, 3L)
      .r4vn_mx_add_directive(store, variable, "factor")
      return(as.name(variable))
    }

    # c.var: continuous.
    if (startsWith(txt, "c.") && nchar(txt) > 2L) {
      variable <- substring(txt, 3L)
      .r4vn_mx_add_directive(store, variable, "continuous")
      return(as.name(variable))
    }

    return(x)
  }

  if (is.call(x)) {
    z <- as.list(x)
    if (length(z) >= 2L) {
      for (j in 2:length(z)) {
        z[[j]] <- .r4vn_mx_translate_expr(z[[j]], store)
      }
    }
    return(as.call(z))
  }

  x
}

.r4vn_mx_is_vars_call <- function(expr) {
  if (!is.call(expr) || !length(expr)) return(FALSE)
  head <- expr[[1L]]

  if (is.symbol(head)) {
    return(identical(as.character(head), "vars"))
  }

  if (is.call(head) &&
      length(head) >= 3L &&
      as.character(head[[1L]]) %in% c("::", ":::")) {
    return(identical(as.character(head[[3L]]), "vars"))
  }

  FALSE
}

.r4vn_mx_unpack_vars <- function(vars_expr, env = parent.frame()) {
  if (is.null(vars_expr) || identical(vars_expr, quote(NULL))) {
    return(list())
  }

  # Inline vars(...) is intentionally not evaluated. This is what allows:
  # vars(c.age, ib2.job*i.treatment)
  # without changing the public table-oriented vars() function.
  if (.r4vn_mx_is_vars_call(vars_expr)) {
    out <- as.list(vars_expr)[-1L]
    if (!length(out)) {
      stop("`vars = vars(...)` must contain at least one model term.", call. = FALSE)
    }
    return(out)
  }

  value <- tryCatch(eval(vars_expr, envir = env), error = function(e) NULL)

  if (inherits(value, "r4vn_vars")) {
    if (!nrow(value)) {
      stop("`vars` contains no model variables.", call. = FALSE)
    }

    specification <- as.character(value$specification)
    bad <- startsWith(specification, "-") |
      grepl("*", specification, fixed = TRUE) |
      specification == "."

    if (any(bad)) {
      stop(
        paste0(
          "Deferred selectors, wildcards, and exclusions are not supported in regression `vars=`. ",
          "Supply explicit model variables."
        ),
        call. = FALSE
      )
    }

    return(lapply(specification, as.name))
  }

  stop(
    "`vars` must be written as `vars(...)` or be an existing `r4vn_vars` object.",
    call. = FALSE
  )
}

.r4vn_mx_build_formula <- function(lhs_expr,
                                   rhs_exprs = list(),
                                   vars_expr = NULL,
                                   env = parent.frame(),
                                   noconstant = FALSE) {
  store <- .r4vn_mx_directive_store()

  lhs_value <- tryCatch(eval(lhs_expr, envir = env), error = function(e) NULL)
  is_formula <- .r4vn_is_formula_expr(lhs_expr) ||
    inherits(lhs_value, "formula")

  vars_items <- .r4vn_mx_unpack_vars(vars_expr, env)

  if (is_formula) {
    if (length(rhs_exprs) || length(vars_items)) {
      stop(
        "When `y` is supplied as a formula, do not also supply predictors through `...` or `vars=`.",
        call. = FALSE
      )
    }

    f <- if (inherits(lhs_value, "formula")) {
      lhs_value
    } else if (inherits(lhs_expr, "formula")) {
      lhs_expr
    } else {
      stats::as.formula(lhs_expr, env = env)
    }

    if (length(f) < 3L) {
      stop("A regression formula must contain both an outcome and predictors.", call. = FALSE)
    }

    f[[3L]] <- .r4vn_mx_translate_expr(f[[3L]], store)

    if (isTRUE(noconstant)) {
      f <- stats::update.formula(f, . ~ . - 1)
    }

    attr(f, "r4vn_model_directives") <- .r4vn_mx_collect_directives(store)
    return(f)
  }

  all_rhs <- c(rhs_exprs, vars_items)
  all_rhs <- .r4vn_mx_drop_redundant_plain(all_rhs)
  lhs <- .r4vn_mx_deparse1(lhs_expr)

  rhs <- if (length(all_rhs)) {
    translated <- lapply(
      all_rhs,
      .r4vn_mx_translate_expr,
      store = store
    )
    paste(
      vapply(translated, .r4vn_mx_deparse1, character(1)),
      collapse = " + "
    )
  } else {
    "1"
  }

  if (isTRUE(noconstant)) {
    rhs <- paste0(rhs, " - 1")
  }

  f <- stats::as.formula(
    paste(lhs, "~", rhs),
    env = env
  )
  attr(f, "r4vn_model_directives") <- .r4vn_mx_collect_directives(store)
  f
}

.r4vn_mx_effective_ref <- function(ref, formula) {
  if (is.null(ref)) return(NULL)

  directives <- attr(formula, "r4vn_model_directives")
  if (is.null(directives) || !length(directives) ||
      !is.list(ref) || is.null(names(ref))) {
    return(ref)
  }

  explicit <- names(directives)[
    vapply(
      directives,
      function(z) !is.null(z$ref),
      logical(1)
    )
  ]

  if (!length(explicit)) return(ref)

  out <- ref[setdiff(names(ref), explicit)]
  if (!length(out)) NULL else out
}


.r4vn_mx_apply_directives <- function(data, formula) {
  directives <- attr(formula, "r4vn_model_directives")
  if (is.null(directives) || !length(directives)) {
    return(data)
  }

  for (nm in names(directives)) {
    d <- directives[[nm]]

    if (!nm %in% names(data)) {
      stop(
        sprintf("Variable `%s` was not found in the model data.", nm),
        call. = FALSE
      )
    }

    x <- data[[nm]]

    if (identical(d$type, "continuous")) {
      if (!is.numeric(x)) {
        stop(
          sprintf(
            "`c.%s` declares `%s` as continuous, but the variable is not numeric.",
            nm, nm
          ),
          call. = FALSE
        )
      }
      next
    }

    f <- if (is.factor(x)) droplevels(x) else factor(x)

    if (nlevels(f) < 2L) {
      stop(
        sprintf(
          "Categorical variable `%s` must have at least two observed levels.",
          nm
        ),
        call. = FALSE
      )
    }

    if (!is.null(d$ref)) {
      lev <- levels(f)

      if (identical(d$ref_mode, "index")) {
        idx <- as.integer(d$ref)
        if (!is.finite(idx) || idx < 1L || idx > length(lev)) {
          stop(
            sprintf(
              "Reference index %s is not available for `%s`. Observed levels: %s.",
              d$ref, nm, paste(lev, collapse = ", ")
            ),
            call. = FALSE
          )
        }
        ref_level <- lev[idx]

      } else if (identical(d$ref_mode, "hybrid")) {
        requested <- as.character(d$ref)[1L]

        if (requested %in% lev) {
          ref_level <- requested
        } else {
          idx <- suppressWarnings(as.integer(requested))

          if (is.finite(idx) && idx >= 1L && idx <= length(lev)) {
            ref_level <- lev[idx]
          } else {
            stop(
              sprintf(
                "Reference `%s` is not available for `%s`. Observed levels: %s.",
                requested, nm, paste(lev, collapse = ", ")
              ),
              call. = FALSE
            )
          }
        }

      } else {
        ref_level <- as.character(d$ref)[1L]
        if (!ref_level %in% lev) {
          stop(
            sprintf(
              "Reference `%s` is not available for `%s`. Observed levels: %s.",
              ref_level, nm, paste(lev, collapse = ", ")
            ),
            call. = FALSE
          )
        }
      }

      f <- stats::relevel(f, ref = ref_level)
    }

    data[[nm]] <- f
  }

  data
}

.r4vn_mx_term_labels <- function(fit) {
  out <- tryCatch(
    attr(stats::terms(fit), "term.labels"),
    error = function(e) character(0)
  )
  as.character(out)
}

.r4vn_mx_clean_term <- function(x) {
  x <- as.character(x)
  x <- gsub("factor\\(([^()]*)\\)", "\\1", x, perl = TRUE)
  x <- gsub(":", " x ", x, fixed = TRUE)
  x
}


# ============================================================================
# Function source: corr.R
# ============================================================================
#' Correlation matrix
#'
#' Computes Pearson, Spearman, or Kendall correlations from variables in a data
#' frame. With `data = NULL`, the active data set is used.
#'
#' @param ... Numeric variables. If omitted, all numeric variables are used.
#' @param data Data frame or `NULL` for active data.
#' @param method Correlation method.
#' @param missing Pairwise or listwise deletion.
#' @param sig Show a p-value matrix.
#' @param obs Show a matrix of pairwise sample sizes.
#' @param ci Show pairwise confidence intervals where available.
#' @param star Add significance stars to the displayed correlation matrix.
#' @param level Confidence level.
#' @param digits,p_digits Decimal places.
#' @param show Logical; open the formatted result in the Viewer. Default `TRUE`.
#' @param console Logical; also print the traditional result in the Console. Default `FALSE`.
#' @return An object of class `r4vn_stat`, returned invisibly. Its `sections`
#'   component contains the formatted correlation matrix and any requested
#'   p-value, pairwise sample-size, or confidence-interval tables. Its `raw`
#'   component contains the numeric correlation (`correlation`), p-value
#'   (`p.value`), and pairwise sample-size (`n`) matrices plus the selected
#'   correlation `method`; `call` records the matched function call.
#' @examples
#'
#' # Extended usage examples
#' d <- data.frame(age = c(20, 25, 30, 35, 40, 45),
#'                 bmi = c(20, 22, 24, 23, 26, 28),
#'                 score = c(60, 65, 68, 72, 75, 80))
#'
#' corr(age, bmi, score, data = d, show = FALSE)
#' corr(age, bmi, score, data = d, method = "spearman", show = FALSE)
#' corr(age, bmi, score, data = d, sig = TRUE, obs = TRUE, ci = TRUE, show = FALSE)
#' corr(age, bmi, score, data = d, star = TRUE, show = FALSE)
#' @export
corr <- function(..., data = NULL, method = c("pearson", "spearman", "kendall"),
                 missing = c("pairwise", "listwise"), sig = FALSE, obs = FALSE,
                 ci = FALSE, star = FALSE, level = 0.95, digits = 3,
                 p_digits = 3, show = TRUE, console = FALSE) {
  method <- match.arg(method); missing <- match.arg(missing)
  if (!is.numeric(level) || length(level) != 1L || level <= 0 || level >= 1) stop("`level` must be between 0 and 1.", call. = FALSE)
  env <- parent.frame(); d <- .r4vn_stat_data(data)
  exprs <- as.list(substitute(list(...)))[-1L]
  if (!length(exprs)) {
    use <- vapply(d, is.numeric, logical(1))
    if (!any(use)) stop("No numeric variables were found.", call. = FALSE)
    vals <- d[use]; nms <- names(vals)
  } else if (length(exprs) == 1L && (.r4vn_is_formula_expr(exprs[[1L]]) || inherits(tryCatch(eval(exprs[[1L]], env), error = function(e) NULL), "formula"))) {
    fv <- tryCatch(eval(exprs[[1L]], env), error = function(e) NULL)
    f <- if (inherits(fv, "formula")) fv else stats::as.formula(exprs[[1L]], env = env)
    vals <- stats::model.frame(f, data = d, na.action = stats::na.pass)
    nms <- names(vals)
  } else {
    vals <- lapply(seq_along(exprs), function(i) .r4vn_eval_var(exprs[[i]], d, env, paste0("variable ", i)))
    nms <- vapply(exprs, .r4vn_deparse1, character(1)); names(vals) <- nms
  }
  if (!all(vapply(vals, is.numeric, logical(1)))) stop("All variables supplied to `corr()` must be numeric.", call. = FALSE)
  X <- as.data.frame(vals, check.names = FALSE)
  if (ncol(X) < 2L) stop("At least two numeric variables are required.", call. = FALSE)
  if (missing == "listwise") X <- X[stats::complete.cases(X), , drop = FALSE]
  k <- ncol(X); R <- P <- N <- matrix(NA_real_, k, k, dimnames = list(names(X), names(X)))
  cil <- list(); idx <- 0L
  for (i in seq_len(k)) for (j in i:k) {
    ok <- stats::complete.cases(X[[i]], X[[j]])
    xi <- X[[i]][ok]; xj <- X[[j]][ok]; nij <- length(xi)
    N[i, j] <- N[j, i] <- nij
    if (i == j) { R[i, j] <- 1; P[i, j] <- 0; next }
    if (nij < 3L || stats::sd(xi) == 0 || stats::sd(xj) == 0) next
    ct <- tryCatch(suppressWarnings(stats::cor.test(xi, xj, method = method, conf.level = level, exact = FALSE)), error = function(e) NULL)
    if (is.null(ct)) next
    r <- unname(ct$estimate); R[i, j] <- R[j, i] <- r; P[i, j] <- P[j, i] <- ct$p.value
    if (ci) {
      idx <- idx + 1L; cint <- if (!is.null(ct$conf.int)) unname(ct$conf.int) else c(NA_real_, NA_real_)
      cil[[idx]] <- data.frame(Variable1 = names(X)[i], Variable2 = names(X)[j], n = nij,
                               Correlation = .r4vn_num(r, digits), Lower = .r4vn_num(cint[1], digits),
                               Upper = .r4vn_num(cint[2], digits), p = .r4vn_p(ct$p.value, p_digits), stringsAsFactors = FALSE)
    }
  }
  disp <- matrix(.r4vn_num(R, digits), nrow = k, dimnames = dimnames(R))
  if (star) {
    marks <- ifelse(P < .001, "***", ifelse(P < .01, "**", ifelse(P < .05, "*", "")))
    marks[is.na(marks)] <- ""; diag(marks) <- ""; disp <- matrix(paste0(.r4vn_num(R, digits), marks), nrow = k, dimnames = dimnames(R))
  }
  sections <- list("Correlation matrix" = disp)
  if (sig) {
    Pshow <- P; diag(Pshow) <- NA_real_
    sections[["P-values"]] <- matrix(.r4vn_p(Pshow, p_digits), nrow = k, dimnames = dimnames(Pshow))
  }
  if (obs) sections[["Observations"]] <- N
  if (ci && length(cil)) sections[["Pairwise confidence intervals"]] <- do.call(rbind, cil)
  note <- if (star) "* p<0.05, ** p<0.01, *** p<0.001." else NULL
  .r4vn_show(.r4vn_result(paste(tools::toTitleCase(method), "correlations"), sections, note,
                           raw = list(correlation = R, p.value = P, n = N, method = method), call = match.call()), show)
}


# ============================================================================
# Function source: regress.R
# ============================================================================
#' Linear regression
#'
#' Fits an ordinary least-squares model. R4VN compact syntax allows models such
#' as `regress(y, c.age, i.sex, i.sex*i.treatment)` without `~` or `+`.
#'
#' @param y Formula or numeric outcome variable.
#' @param ... Predictors or model terms when `y` is not a formula.
#' @param vars Optional model terms written as `vars(...)`.
#' @param data Data frame or `NULL` for active data.
#' @param noconstant Fit without an intercept.
#' @param vce Model-based, HC1 robust, or cluster-robust covariance.
#' @param cluster Cluster variable used when `vce = "cluster"`.
#' @param weights Optional non-negative analytic weights.
#' @param subset Optional logical subset expression.
#' @param ref Optional named list of factor reference levels.
#' @param standardized Also display standardized coefficients for numeric columns.
#' @param vif Also display coefficient-level variance inflation factors.
#' @param diagnosis Logical; if `TRUE`, append model-diagnostic tables. For linear regression these include residual normality, a Breusch-Pagan heteroscedasticity test, standardized/studentized residuals, leverage, Cook's distance, DFFITS, COVRATIO, influential observations, and collinearity diagnostics. Default `FALSE`.
#' @param level Confidence level.
#' @param digits,p_digits Decimal places.
#' @param show Logical; open the formatted result in the Viewer. Default `TRUE`.
#' @param console Logical; also print the traditional result in the Console. Default `FALSE`.
#'
#' @return An object of class `r4vn_stat`, returned invisibly. Its `sections`
#'   component contains the formatted model summary, ANOVA, coefficient table,
#'   and any requested standardized-coefficient or VIF tables. In `raw`,
#'   `model` is the fitted `lm` object, `vcov` is the covariance matrix,
#'   `coefficients` contains coefficient-level estimates and tests, `overall`
#'   contains the overall model test, and `vce` and `model.terms` record the
#'   covariance estimator and fitted terms.
#'
#' @details
#' Compact model prefixes are `c.x` for continuous, `i.x` for categorical,
#' `b2.x` for the second factor level as reference, and `ib2.x` for value/level
#' 2 as reference. Use `*` for main effects plus interaction and `:` for
#' interaction only.
#'
#' Ordinary linear regression models should be compared with the usual nested
#' F test rather than [lrtest()].
#'
#' @examples
#' d <- data.frame(score = c(60, 65, 68, 72, 75, 80, 77, 70),
#'                 age = c(20, 25, 30, 35, 40, 45, 50, 55),
#'                 bmi = c(20, 22, 24, 23, 26, 28, 27, 25),
#'                 sex = factor(rep(c("Female", "Male"), 4)))
#'
#' regress(score, age, bmi, i.sex, data = d, show = FALSE)
#' regress(score, vars = vars(c.age, c.bmi, i.sex), data = d, show = FALSE)
#' # Full model diagnostics
#' m <- regress(score, c.age, c.bmi, i.sex, data = d, diagnosis = TRUE, show = FALSE)
#' m$sections$`Model diagnosis`
#' m$sections$`Influence diagnostics`
#' # Postestimation diagnostics can also be generated as variables
#' usedf(d)
#' regress(score, c.age, c.bmi, i.sex, diagnosis = FALSE, show = FALSE)
#' predict(newvar = stdres, type = "standardized", show = FALSE)
#' predict(newvar = cooksd, type = "cooksd", show = FALSE)
#' @export
regress <- function(y, ..., vars = NULL, data = NULL, noconstant = FALSE,
                    vce = c("model", "robust", "cluster"), cluster = NULL,
                    weights = NULL, subset = NULL, ref = NULL,
                    standardized = FALSE, vif = FALSE, diagnosis = FALSE, level = 0.95,
                    digits = 3, p_digits = 3, show = TRUE, console = FALSE) {
  vce <- match.arg(vce); env <- parent.frame()
  if (!is.numeric(level) || length(level) != 1L || level <= 0 || level >= 1) stop("`level` must be between 0 and 1.", call. = FALSE)

  rhs <- as.list(substitute(list(...)))[-1L]
  vars_expr <- if (missing(vars)) NULL else substitute(vars)
  f <- .r4vn_mx_build_formula(substitute(y), rhs, vars_expr, env, noconstant)

  prep <- .r4vn_prepare_model_data(
    data, env,
    substitute(subset),
    substitute(weights),
    substitute(cluster),
    ref = .r4vn_mx_effective_ref(ref, f)
  )
  prep$data <- .r4vn_mx_apply_directives(prep$data, f)

  fit_formula <- .r4vn_model_formula(f, nrow(prep$data), names(prep$data), weights = prep$weights)
  fit <- stats::lm(fit_formula, data = prep$data, weights = .r4vn_internal_weights_7e4f9c,
                   na.action = stats::na.omit, model = TRUE, x = TRUE, y = TRUE)
  fit <- .r4vn_mx_add_reference_aliases(fit, f)
  used <- .r4vn_used_rows(fit, nrow(prep$data)); clu <- if (is.null(prep$cluster)) NULL else prep$cluster[used]
  V <- .r4vn_model_vcov(fit, vce, clu)
  cr <- .r4vn_coef_raw(fit, V, level, "t")
  sm <- summary(fit); a <- stats::anova(fit); n <- stats::nobs(fit)
  if (vce == "model" && !is.null(sm$fstatistic)) {
    overall <- list(statistic = unname(sm$fstatistic[1]), df1 = unname(sm$fstatistic[2]), df2 = unname(sm$fstatistic[3]),
                    p.value = stats::pf(sm$fstatistic[1], sm$fstatistic[2], sm$fstatistic[3], lower.tail = FALSE))
  } else overall <- .r4vn_wald_overall(fit, V, linear = TRUE)
  dep <- .r4vn_deparse1(f[[2L]])
  info <- data.frame(Statistic = c("Dependent variable", "Number of obs", sprintf("F(%s, %s)", overall$df1, overall$df2), "Prob > F", "R-squared", "Adjusted R-squared", "Root MSE", "VCE"),
                     Value = c(dep, n, .r4vn_num(overall$statistic, 2), .r4vn_p(overall$p.value, p_digits),
                               .r4vn_num(sm$r.squared, digits), .r4vn_num(sm$adj.r.squared, digits),
                               .r4vn_num(sm$sigma, digits), vce), stringsAsFactors = FALSE)
  acol <- function(nm) if (nm %in% names(a)) a[[nm]] else rep(NA_real_, nrow(a))
  av <- data.frame(Source = rownames(a), SS = .r4vn_num(acol("Sum Sq"), digits), df = acol("Df"),
                   MS = .r4vn_num(acol("Mean Sq"), digits), F = .r4vn_num(acol("F value"), 2),
                   p = .r4vn_p(acol("Pr(>F)"), p_digits), stringsAsFactors = FALSE, check.names = FALSE)
  sections <- list("Model summary" = info, "ANOVA" = av,
                   "Coefficients" = .r4vn_coef_table(cr, digits, p_digits, "t", FALSE, "Coefficient"))
  if (standardized) {
    mf <- stats::model.frame(fit); yy <- stats::model.response(mf); X <- stats::model.matrix(fit)
    sx <- apply(X, 2, stats::sd); sy <- stats::sd(yy); beta <- stats::coef(fit) * sx / sy; beta[names(beta) == "(Intercept)"] <- NA_real_
    sections[["Standardized coefficients"]] <- data.frame(Term = sub("^\\(Intercept\\)$", "_cons", names(beta)), Beta = .r4vn_num(beta, digits), stringsAsFactors = FALSE)
  }
  if (vif) { vv <- .r4vn_vif(fit, digits); if (!is.null(vv)) sections[["Variance inflation factors"]] <- vv }
  diagnostics <- NULL
  if (isTRUE(diagnosis)) {
    diagnostics <- .r4vn_model_diagnosis(fit, kind = "linear", digits = digits, p_digits = p_digits)
    sections <- c(sections, diagnostics)
  }
  .r4vn_show(.r4vn_result("Linear regression", sections,
                           if (vce == "model") NULL else "The model F test uses the requested robust covariance; the SS/MS table remains the ordinary least-squares decomposition.",
                           raw = list(model = fit, vcov = V, coefficients = cr, overall = overall, vce = vce,
                                      model.terms = .r4vn_mx_term_labels(fit), diagnostics = diagnostics), call = match.call()), show)
}


# ============================================================================
# Function source: logistic.R
# ============================================================================
#' Binary logistic regression
#'
#' Fits binary logistic regression using formula syntax or compact R4VN syntax.
#' Compact syntax avoids the need to type `~` and `+`.
#'
#' @param y Formula or binary outcome variable.
#' @param ... Predictors or model terms when `y` is not a formula.
#' @param vars Optional model terms written as `vars(...)`. The expression is
#'   captured without evaluating the public table-oriented `vars()` parser, so
#'   interactions are allowed.
#' @param data Data frame or `NULL` for active data.
#' @param event Event level for a simple named outcome.
#' @param or Add an odds-ratio table while retaining coefficients.
#' @param exp Display odds ratios only.
#' @param noconstant Fit without an intercept.
#' @param vce Model-based, HC1 robust, or cluster-robust covariance.
#' @param cluster Cluster variable.
#' @param weights Optional non-negative weights.
#' @param subset Optional logical subset.
#' @param ref Optional named list of factor reference levels.
#' @param gof Show a Hosmer-Lemeshow test.
#' @param groups Number of groups for the Hosmer-Lemeshow test.
#' @param classification Show a classification table.
#' @param cutoff Classification cutoff.
#' @param vif Show coefficient-level VIFs.
#' @param diagnosis Logical; if `TRUE`, append calibration/goodness-of-fit, discrimination, residual, influence, and collinearity diagnostics appropriate for binary logistic regression. Default `FALSE`.
#' @param level Confidence level.
#' @param digits,p_digits Decimal places.
#' @param show Logical; open the formatted result in the Viewer. Default `TRUE`.
#' @param console Logical; also print the traditional result in the Console. Default `FALSE`.
#'
#' @return An object of class `r4vn_stat`, returned invisibly. Its `sections`
#'   component contains the formatted model summary, coefficient and/or odds-
#'   ratio tables, and any requested goodness-of-fit, classification, or VIF
#'   tables. In `raw`, `model` is the fitted binomial `glm` object, `vcov` is
#'   the covariance matrix, `coefficients` contains coefficient-level estimates
#'   and tests, `logLik` and `null.logLik` are model log likelihoods, `pseudo.r2`
#'   is McFadden-style pseudo-R-squared, `event` records the modeled outcome
#'   level, and `vce` and `model.terms` record the covariance estimator and
#'   fitted terms.
#'
#' @details
#' Compact model syntax:
#'
#' \itemize{
#'   \item `x`: use the variable as stored in the data.
#'   \item `c.x`: force `x` to be continuous.
#'   \item `i.x`: force `x` to be categorical.
#'   \item `b2.x`, `b3.x`, ...: categorical with the corresponding factor-level
#'         position as reference.
#'   \item `ib0.x`, `ib1.x`, `ib2.x`, ...: categorical with the requested
#'         value/level as reference. If that literal level is unavailable, a
#'         positive integer can fall back to the corresponding factor-level
#'         position.
#'   \item `i.a*i.b`: main effects for `a` and `b` plus their interaction.
#'   \item `i.a:i.b`: interaction only.
#'   \item `c.x*i.a`: continuous and categorical main effects plus interaction.
#' }
#'
#' Thus `logistic(y, ib2.occupation*i.treatment, c.age)` fits occupation,
#' treatment, occupation-by-treatment interaction, and age without requiring
#' formula operators `~` or `+`.
#'
#' The fitted `glm` object is stored in `result$raw$model`, so nested models can
#' be compared directly with [lrtest()].
#'
#' @examples
#' set.seed(2026)
#' d <- data.frame(
#'   outcome = factor(rbinom(200, 1, .35), levels = 0:1,
#'                    labels = c("No", "Yes")),
#'   age = rnorm(200, 45, 12),
#'   occupation = factor(sample(c("Office", "Worker", "Other"), 200, TRUE)),
#'   treatment = factor(sample(c("No", "Yes"), 200, TRUE))
#' )
#'
#' m1 <- logistic(
#'   outcome,
#'   c.age,
#'   i.occupation,
#'   i.treatment,
#'   data = d,
#'   event = "Yes",
#'   show = FALSE
#' )
#'
#' m2 <- logistic(
#'   outcome,
#'   c.age,
#'   ib2.occupation*i.treatment,
#'   data = d,
#'   event = "Yes",
#'   show = FALSE
#' )
#'
#' m3 <- logistic(
#'   outcome,
#'   vars = vars(c.age, ib2.occupation*i.treatment),
#'   data = d,
#'   event = "Yes",
#'   show = FALSE
#' )
#'
#' lrtest(m1, m2, show = FALSE)
#'
#' # Request a complete diagnostic panel
#' logistic(outcome, c.age, i.occupation, data = d, event = "Yes",
#'          diagnosis = TRUE, show = FALSE)
#'
#' @seealso [lrtest()], [poisson()]
#' @export
logistic <- function(y, ..., vars = NULL, data = NULL, event = NULL,
                     or = FALSE, exp = FALSE,
                     noconstant = FALSE, vce = c("model", "robust", "cluster"),
                     cluster = NULL, weights = NULL, subset = NULL, ref = NULL,
                     gof = FALSE, groups = 10, classification = FALSE, cutoff = 0.5,
                     vif = FALSE, diagnosis = FALSE, level = 0.95, digits = 3, p_digits = 3,
                     show = TRUE, console = FALSE) {
  vce <- match.arg(vce); env <- parent.frame()
  rhs <- as.list(substitute(list(...)))[-1L]
  vars_expr <- if (missing(vars)) NULL else substitute(vars)

  if (!is.numeric(level) || length(level) != 1L || level <= 0 || level >= 1) stop("`level` must be between 0 and 1.", call. = FALSE)
  if (!is.numeric(cutoff) || length(cutoff) != 1L || cutoff <= 0 || cutoff >= 1) stop("`cutoff` must be between 0 and 1.", call. = FALSE)
  if (!is.numeric(groups) || length(groups) != 1L || groups < 3) stop("`groups` must be at least 3.", call. = FALSE)

  f <- .r4vn_mx_build_formula(substitute(y), rhs, vars_expr, env, noconstant)

  prep <- .r4vn_prepare_model_data(
    data, env,
    substitute(subset),
    substitute(weights),
    substitute(cluster),
    ref = .r4vn_mx_effective_ref(ref, f)
  )
  prep$data <- .r4vn_mx_apply_directives(prep$data, f)

  event_label <- NULL; lhs <- f[[2L]]
  if (is.symbol(lhs) && as.character(lhs) %in% names(prep$data)) {
    nm <- as.character(lhs); vv <- prep$data[[nm]]; lev <- unique(as.character(vv[!is.na(vv)]))
    if (length(lev) != 2L) stop("The logistic outcome must have exactly two observed values.", call. = FALSE)
    if (!is.null(event)) {
      event_label <- as.character(event)[1L]
      if (!event_label %in% lev) stop("`event` was not found in the outcome.", call. = FALSE)
      prep$data[[nm]] <- as.integer(as.character(vv) == event_label)
    } else if (is.factor(vv) || is.character(vv)) {
      prep$data[[nm]] <- droplevels(factor(vv))
      event_label <- levels(prep$data[[nm]])[2L]
    } else if (is.logical(vv)) {
      event_label <- "TRUE"
      prep$data[[nm]] <- as.integer(vv)
    } else if (is.numeric(vv) && all(lev %in% c("0", "1"))) {
      event_label <- "1"
    } else stop("A numeric logistic outcome must be coded 0 and 1, or specify `event`.", call. = FALSE)
  } else if (!is.null(event)) stop("`event` can be used only when the formula outcome is a simple variable name.", call. = FALSE)

  fit_formula <- .r4vn_model_formula(f, nrow(prep$data), names(prep$data), weights = prep$weights)
  fit <- stats::glm(fit_formula, data = prep$data, family = stats::binomial(),
                    weights = .r4vn_internal_weights_7e4f9c,
                    na.action = stats::na.omit, model = TRUE, x = TRUE, y = TRUE)
  fit <- .r4vn_mx_add_reference_aliases(fit, f)
  used <- .r4vn_used_rows(fit, nrow(prep$data)); clu <- if (is.null(prep$cluster)) NULL else prep$cluster[used]
  V <- .r4vn_model_vcov(fit, vce, clu); cr <- .r4vn_coef_raw(fit, V, level, "z")
  nullfit <- .r4vn_glm_null(fit)
  ll <- as.numeric(stats::logLik(fit)); ll0 <- if (is.null(nullfit)) NA_real_ else as.numeric(stats::logLik(nullfit))
  if (vce == "model") {
    chi <- fit$null.deviance - fit$deviance; df <- fit$df.null - fit$df.residual; pp <- stats::pchisq(chi, df, lower.tail = FALSE); test_name <- "LR chi2"
  } else {
    ov <- .r4vn_wald_overall(fit, V, FALSE); chi <- ov$statistic; df <- ov$df1; pp <- ov$p.value; test_name <- "Wald chi2"
  }
  pseudo <- if (is.finite(ll0) && ll0 != 0) 1 - ll / ll0 else NA_real_
  dep <- .r4vn_deparse1(f[[2L]])
  info <- data.frame(Statistic = c("Dependent variable", "Number of obs", "Events", sprintf("%s(%s)", test_name, df), "Prob > chi2", "Log likelihood", "Pseudo R2", "AIC", "BIC", "VCE"),
                     Value = c(dep, stats::nobs(fit), sum(fit$y), .r4vn_num(chi, 2), .r4vn_p(pp, p_digits),
                               .r4vn_num(ll, digits), .r4vn_num(pseudo, digits), .r4vn_num(stats::AIC(fit), digits),
                               .r4vn_num(stats::BIC(fit), digits), vce), stringsAsFactors = FALSE)
  sections <- list("Model summary" = info)
  if (!isTRUE(exp)) sections[["Coefficients"]] <- .r4vn_coef_table(cr, digits, p_digits, "z", FALSE, "Coefficient")
  if (isTRUE(exp) || isTRUE(or)) sections[["Odds ratios"]] <- .r4vn_coef_table(cr, digits, p_digits, "z", TRUE, "Odds ratio")
  if (gof) sections[["Goodness of fit"]] <- .r4vn_logistic_gof(fit, as.integer(groups), digits, p_digits)
  if (classification) {
    cc <- .r4vn_classification(fit, cutoff, 1); sections[["Classification table"]] <- cc$table; sections[["Classification statistics"]] <- cc$metrics
  }
  if (vif) { vv <- .r4vn_vif(fit, digits); if (!is.null(vv)) sections[["Variance inflation factors"]] <- vv }
  diagnostics <- NULL
  if (isTRUE(diagnosis)) {
    diagnostics <- .r4vn_model_diagnosis(fit, kind = "logistic", groups = as.integer(groups), digits = digits, p_digits = p_digits)
    sections <- c(sections, diagnostics)
  }
  note <- if (!is.null(event_label)) paste0("Modeled event: ", event_label, ".") else NULL
  .r4vn_show(.r4vn_result("Logistic regression", sections, note,
                           raw = list(model = fit, vcov = V, coefficients = cr, logLik = ll, null.logLik = ll0,
                                      pseudo.r2 = pseudo, event = event_label, vce = vce,
                                      model.terms = .r4vn_mx_term_labels(fit), diagnostics = diagnostics), call = match.call()), show)
}



# ============================================================================
# Internal Poisson-family compatibility helper
# ============================================================================
.r4vn_poisson_family <- function(link = "log") {
  # Ordinary user-facing character links.
  if (is.character(link) && length(link) == 1L && !is.na(link)) {
    link_name <- match.arg(
      tolower(link),
      c("log", "identity", "sqrt")
    )

    return(do.call(
      stats::poisson,
      list(link = link_name)
    ))
  }

  # MASS::glm.nb() can call an attached/exported poisson() as
  # poisson(link = log), i.e. with the actual function rather than "log".
  if (is.function(link)) {
    if (identical(link, base::log)) {
      return(stats::poisson(link = "log"))
    }
    if (identical(link, base::identity)) {
      return(stats::poisson(link = "identity"))
    }
    if (identical(link, base::sqrt)) {
      return(stats::poisson(link = "sqrt"))
    }
  }

  # Also accept an already constructed link-glm object.
  if (inherits(link, "link-glm")) {
    nm <- tryCatch(link$name, error = function(e) NULL)

    if (is.character(nm) && length(nm) == 1L && !is.na(nm)) {
      nm <- match.arg(
        tolower(nm),
        c("log", "identity", "sqrt")
      )

      return(do.call(
        stats::poisson,
        list(link = nm)
      ))
    }
  }

  stop(
    "`link` must be one of \"log\", \"identity\", or \"sqrt\", ",
    "or the corresponding link function.",
    call. = FALSE
  )
}

# ============================================================================
# Function source: poisson.R
# ============================================================================
#' Poisson regression or Poisson family
#'
#' Fits Poisson regression using formula syntax or compact R4VN syntax. When
#' called without a model, for example `poisson()` or
#' `poisson(link = "identity")`, returns the ordinary [stats::poisson()] family.
#' A binary outcome is also supported; `event` explicitly identifies the event
#' category and `rr = TRUE` requests a risk-ratio display. Robust or
#' cluster-robust VCE is generally appropriate for modified-Poisson binary models.
#'
#' @usage
#' poisson(
#'   y, ..., vars = NULL, data = NULL, exposure = NULL, offset = NULL,
#'   event = NULL, irr = FALSE, rr = FALSE, exp = FALSE, link = "log",
#'   noconstant = FALSE, vce = c("model", "robust", "cluster"), cluster = NULL,
#'   weights = NULL, subset = NULL, ref = NULL, vif = FALSE, diagnosis = FALSE,
#'   level = 0.95, digits = 3, p_digits = 3, show = TRUE, console = FALSE
#' )
#'
#' @param y Formula or count outcome. Omit to obtain the base R Poisson family.
#' @param ... Predictors or model terms when `y` is not a formula.
#' @param vars Optional model terms written as `vars(...)`.
#' @param data Data frame or `NULL` for active data.
#' @param exposure Optional person-time variable; its logarithm is used as offset.
#' @param offset Optional offset already on the linear-predictor scale.
#' @param event Event value when `y` is binary. For a 0/1 variable the default
#'   event is 1; for a factor, the second level is used unless specified.
#' @param irr Add an incidence-rate-ratio table for count outcomes.
#' @param rr Add a risk-ratio table for binary outcomes.
#' @param exp Display exponentiated coefficients only (IRR for counts, RR for binary outcomes).
#' @param link Link used only in family mode. Accepts `"log"`, `"identity"`,
#'   or `"sqrt"`; the corresponding link functions are also accepted for
#'   compatibility with packages such as MASS.
#' @param noconstant Fit without an intercept.
#' @param vce Model-based, HC1 robust, or cluster-robust covariance.
#' @param cluster Cluster variable.
#' @param weights Optional non-negative weights.
#' @param subset Optional logical subset.
#' @param ref Optional named list of factor reference levels.
#' @param vif Show coefficient-level VIFs.
#' @param diagnosis Logical; if `TRUE`, append Poisson model diagnostics including Pearson/deviance dispersion, goodness-of-fit, residual/influence measures, influential observations, and collinearity diagnostics. Default `FALSE`.
#' @param level Confidence level.
#' @param digits,p_digits Decimal places.
#' @param show Logical; open the formatted result in the Viewer. Default `TRUE`.
#' @param console Logical; also print the traditional result in the Console. Default `FALSE`.
#'
#' @return If `y` is omitted, a base-R `family` object for the Poisson
#'   distribution with the requested link is returned for compatibility with
#'   modeling functions. Otherwise an object of class `r4vn_stat` is returned
#'   invisibly. Its `sections` component contains the formatted model summary,
#'   coefficient and/or exponentiated-effect tables, goodness-of-fit results,
#'   and any requested VIF table. In `raw`, `model` is the fitted Poisson `glm`
#'   object, `vcov` is the covariance matrix, `coefficients` contains
#'   coefficient-level estimates and tests, `logLik` and `null.logLik` are model
#'   log likelihoods, `pearson` is the Pearson chi-square statistic, `offset`
#'   stores the offset used by the fitted model, and `event`, `binary`, `vce`,
#'   and `model.terms` describe binary-event handling, covariance estimation,
#'   and fitted terms.
#'
#' @details
#' The same compact syntax used by [logistic()] is supported: `c.x`, `i.x`,
#' `b2.x`, `ib2.x`, `*` for main effects plus interaction, and `:` for
#' interaction only.
#'
#' Standard Poisson models fitted with model-based VCE can be compared using
#' [lrtest()]. Quasi-Poisson models do not have an ordinary likelihood and are
#' not supported by `lrtest()`.
#'
#' @examples
#' set.seed(2026)
#' d <- data.frame(
#'   cases = rpois(200, 2),
#'   time = runif(200, .5, 4),
#'   age = rnorm(200, 45, 12),
#'   sex = factor(sample(c("Female", "Male"), 200, TRUE)),
#'   treatment = factor(sample(c("No", "Yes"), 200, TRUE))
#' )
#'
#' p1 <- poisson(
#'   cases,
#'   c.age,
#'   i.sex,
#'   i.treatment,
#'   data = d,
#'   exposure = time,
#'   show = FALSE
#' )
#'
#' p2 <- poisson(
#'   cases,
#'   c.age,
#'   i.sex*i.treatment,
#'   data = d,
#'   exposure = time,
#'   show = FALSE
#' )
#'
#' lrtest(p1, p2, show = FALSE)
#'
#' # Modified Poisson for a binary outcome
#' d$event01 <- as.integer(d$cases > 1)
#' poisson(event01, c.age, i.sex, data = d, event = 1, rr = TRUE, vce = "robust")
#'
#' # Dispersion, residual, influence, and collinearity diagnostics
#' poisson(cases, c.age, i.sex, data = d, exposure = time, diagnosis = TRUE, show = FALSE)
#'
#' @seealso [logistic()], [lrtest()]
#' @export
poisson <- function(y, ..., vars = NULL, data = NULL, exposure = NULL, offset = NULL, event = NULL,
                    irr = FALSE, rr = FALSE, exp = FALSE, link = "log", noconstant = FALSE,
                    vce = c("model", "robust", "cluster"), cluster = NULL,
                    weights = NULL, subset = NULL, ref = NULL, vif = FALSE, diagnosis = FALSE,
                    level = 0.95, digits = 3, p_digits = 3,
                    show = TRUE, console = FALSE) {
  if (missing(y)) return(.r4vn_poisson_family(link))
  vce <- match.arg(vce); env <- parent.frame()
  rhs <- as.list(substitute(list(...)))[-1L]
  vars_expr <- if (missing(vars)) NULL else substitute(vars)

  if (!is.numeric(level) || length(level) != 1L || level <= 0 || level >= 1) stop("`level` must be between 0 and 1.", call. = FALSE)
  fam <- .r4vn_poisson_family(link)
  if ((isTRUE(irr) || isTRUE(rr) || isTRUE(exp)) && fam$link != "log") stop("Exponentiated IRR/RR display requires the log link.", call. = FALSE)

  f <- .r4vn_mx_build_formula(substitute(y), rhs, vars_expr, env, noconstant)

  prep <- .r4vn_prepare_model_data(
    data, env,
    substitute(subset),
    substitute(weights),
    substitute(cluster),
    substitute(exposure),
    substitute(offset),
    .r4vn_mx_effective_ref(ref, f)
  )
  prep$data <- .r4vn_mx_apply_directives(prep$data, f)

  # Binary outcomes are valid in a log-link Poisson model (commonly with robust
  # standard errors to estimate risk ratios). `event` makes the modeled event
  # explicit and mirrors logistic(). Count outcomes continue unchanged.
  event_label <- NULL
  binary_outcome <- FALSE
  lhs <- f[[2L]]
  if (is.symbol(lhs) && as.character(lhs) %in% names(prep$data)) {
    nm <- as.character(lhs); vv <- prep$data[[nm]]; lev <- unique(as.character(vv[!is.na(vv)]))
    if (!is.null(event)) {
      if (length(lev) != 2L) stop("`event` can be used only with a binary Poisson outcome.", call. = FALSE)
      event_label <- as.character(event)[1L]
      if (!event_label %in% lev) stop("`event` was not found in the Poisson outcome.", call. = FALSE)
      prep$data[[nm]] <- as.integer(as.character(vv) == event_label)
      binary_outcome <- TRUE
    } else if ((is.factor(vv) || is.character(vv) || is.logical(vv)) && length(lev) == 2L) {
      event_label <- if (is.factor(vv)) levels(droplevels(vv))[2L] else tail(lev, 1L)
      prep$data[[nm]] <- as.integer(as.character(vv) == event_label)
      binary_outcome <- TRUE
    } else if (is.numeric(vv) && length(lev) == 2L && all(lev %in% c("0", "1"))) {
      event_label <- "1"
      binary_outcome <- TRUE
    }
  } else if (!is.null(event)) stop("`event` can be used only when the outcome is a simple variable name.", call. = FALSE)

  if (!is.null(prep$exposure) && !is.null(prep$offset)) stop("Use either `exposure` or `offset`, not both.", call. = FALSE)
  off <- if (!is.null(prep$exposure)) base::log(prep$exposure) else prep$offset
  fit_formula <- .r4vn_model_formula(f, nrow(prep$data), names(prep$data),
                                      weights = prep$weights, offset = off)
  fit <- stats::glm(fit_formula, data = prep$data, family = fam,
                    weights = .r4vn_internal_weights_7e4f9c,
                    offset = .r4vn_internal_offset_7e4f9c,
                    na.action = stats::na.omit, model = TRUE, x = TRUE, y = TRUE)
  fit <- .r4vn_mx_add_reference_aliases(fit, f)
  if (any(fit$y < 0) || any(abs(fit$y - round(fit$y)) > sqrt(.Machine$double.eps))) stop("The Poisson outcome must contain non-negative integer counts.", call. = FALSE)
  used <- .r4vn_used_rows(fit, nrow(prep$data)); clu <- if (is.null(prep$cluster)) NULL else prep$cluster[used]
  V <- .r4vn_model_vcov(fit, vce, clu); cr <- .r4vn_coef_raw(fit, V, level, "z")
  off_used <- if (is.null(off)) NULL else off[used]
  nullfit <- .r4vn_glm_null(fit)
  ll <- as.numeric(stats::logLik(fit)); ll0 <- if (is.null(nullfit)) NA_real_ else as.numeric(stats::logLik(nullfit))
  if (vce == "model") {
    chi <- fit$null.deviance - fit$deviance; df <- fit$df.null - fit$df.residual; pp <- stats::pchisq(chi, df, lower.tail = FALSE); test_name <- "LR chi2"
  } else {
    ov <- .r4vn_wald_overall(fit, V, FALSE); chi <- ov$statistic; df <- ov$df1; pp <- ov$p.value; test_name <- "Wald chi2"
  }
  pearson <- sum(stats::residuals(fit, type = "pearson")^2)
  dep <- .r4vn_deparse1(f[[2L]])
  info <- data.frame(Statistic = c("Dependent variable", "Number of obs", sprintf("%s(%s)", test_name, df), "Prob > chi2", "Log likelihood", "Pseudo R2", "AIC", "BIC", "Deviance", "Pearson chi2", "Pearson dispersion", "VCE"),
                     Value = c(dep, stats::nobs(fit), .r4vn_num(chi, 2), .r4vn_p(pp, p_digits), .r4vn_num(ll, digits),
                               .r4vn_num(if (is.finite(ll0) && ll0 != 0) 1 - ll / ll0 else NA_real_, digits),
                               .r4vn_num(stats::AIC(fit), digits), .r4vn_num(stats::BIC(fit), digits),
                               .r4vn_num(fit$deviance, digits), .r4vn_num(pearson, digits),
                               .r4vn_num(pearson / fit$df.residual, digits), vce), stringsAsFactors = FALSE)
  gof <- data.frame(Test = c("Deviance", "Pearson"), Chi.square = .r4vn_num(c(fit$deviance, pearson), digits),
                    df = fit$df.residual, p = .r4vn_p(stats::pchisq(c(fit$deviance, pearson), fit$df.residual, lower.tail = FALSE), p_digits), stringsAsFactors = FALSE)
  sections <- list("Model summary" = info)
  if (!isTRUE(exp)) sections[["Coefficients"]] <- .r4vn_coef_table(cr, digits, p_digits, "z", FALSE, "Coefficient")
  if (isTRUE(exp) || isTRUE(irr) || isTRUE(rr)) {
    if (isTRUE(binary_outcome)) sections[["Risk ratios"]] <- .r4vn_coef_table(cr, digits, p_digits, "z", TRUE, "Risk ratio")
    else sections[["Incidence-rate ratios"]] <- .r4vn_coef_table(cr, digits, p_digits, "z", TRUE, "IRR")
  }
  sections[["Goodness of fit"]] <- gof
  if (vif) { vv <- .r4vn_vif(fit, digits); if (!is.null(vv)) sections[["Variance inflation factors"]] <- vv }
  diagnostics <- NULL
  if (isTRUE(diagnosis)) {
    diagnostics <- .r4vn_model_diagnosis(fit, kind = "poisson", digits = digits, p_digits = p_digits)
    sections <- c(sections, diagnostics)
  }
  note <- c(
    if (!is.null(prep$exposure)) "The logarithm of exposure was included as an offset." else if (!is.null(prep$offset)) "An offset was included on the linear-predictor scale." else NULL,
    if (isTRUE(binary_outcome)) paste0("Modeled binary event: ", event_label, ". With the log link, exponentiated coefficients are risk ratios; robust/cluster VCE is generally preferred for binary-outcome modified Poisson inference.") else NULL
  )
  .r4vn_show(.r4vn_result("Poisson regression", sections, note,
                           raw = list(model = fit, vcov = V, coefficients = cr, logLik = ll, null.logLik = ll0,
                                      pearson = pearson, offset = off_used, event = event_label, binary = binary_outcome, vce = vce,
                                      model.terms = .r4vn_mx_term_labels(fit), diagnostics = diagnostics), call = match.call()), show)
}


# ============================================================================
# Likelihood-ratio test helpers
# ============================================================================

.r4vn_lr_extract <- function(x, label) {
  vce <- NULL

  if (is.list(x) &&
      !is.null(x$raw) &&
      is.list(x$raw) &&
      !is.null(x$raw$model)) {
    if (!is.null(x$raw$vce)) {
      vce <- as.character(x$raw$vce)[1L]
    }
    x <- x$raw$model
  }

  if (!is.null(vce) && !identical(vce, "model")) {
    stop(
      sprintf(
        paste0(
          "`%s` was fitted with vce = \"%s\". ",
          "Likelihood-ratio comparison requires model-based likelihood inference. ",
          "Refit the compared R4VN models with vce = \"model\"."
        ),
        label, vce
      ),
      call. = FALSE
    )
  }

  if (inherits(x, "lm") && !inherits(x, "glm") &&
      !inherits(x, "survreg")) {
    stop(
      paste0(
        "`lrtest()` is not used for ordinary linear regression in R4VN. ",
        "Use the nested-model F test instead."
      ),
      call. = FALSE
    )
  }

  supported <- inherits(x, "glm") ||
    inherits(x, "negbin") ||
    inherits(x, "coxph") ||
    inherits(x, "survreg")

  if (!supported) {
    stop(
      sprintf(
        "`%s` is not a supported likelihood-based regression model.",
        label
      ),
      call. = FALSE
    )
  }

  list(model = x, vce = vce)
}

.r4vn_lr_signature <- function(model, label) {
  if (inherits(model, "negbin")) {
    link <- tryCatch(model$family$link, error = function(e) "")
    return(list(
      type = "Negative binomial regression",
      signature = paste0("negbin:", link)
    ))
  }

  if (inherits(model, "coxph")) {
    if (!is.null(model$naive.var)) {
      stop(
        sprintf(
          "`%s` is a Cox model using robust variance; a standard likelihood-ratio chi-square comparison is not appropriate.",
          label
        ),
        call. = FALSE
      )
    }

    method <- if (!is.null(model$method)) as.character(model$method)[1L] else ""
    return(list(
      type = "Cox regression",
      signature = paste0("coxph:", method)
    ))
  }

  if (inherits(model, "survreg")) {
    dist <- if (!is.null(model$dist)) as.character(model$dist)[1L] else ""
    return(list(
      type = "Parametric survival regression",
      signature = paste0("survreg:", dist)
    ))
  }

  if (inherits(model, "glm")) {
    fam <- stats::family(model)
    fam_name <- tolower(fam$family)
    link <- fam$link

    if (grepl("^quasi", fam_name)) {
      stop(
        sprintf(
          "`%s` uses a quasi-likelihood family. A classical likelihood-ratio test is not available.",
          label
        ),
        call. = FALSE
      )
    }

    type <- if (identical(fam_name, "binomial") && identical(link, "logit")) {
      "Logistic regression"
    } else if (identical(fam_name, "poisson")) {
      "Poisson regression"
    } else {
      paste0("GLM (", fam$family, ", ", link, " link)")
    }

    return(list(
      type = type,
      signature = paste("glm", fam_name, link, sep = ":")
    ))
  }

  stop(sprintf("`%s` is not a supported model.", label), call. = FALSE)
}

.r4vn_lr_model_frame <- function(model) {
  tryCatch(stats::model.frame(model), error = function(e) NULL)
}

.r4vn_lr_model_matrix <- function(model) {
  tryCatch(stats::model.matrix(model), error = function(e) NULL)
}

.r4vn_lr_model_info <- function(model, label, r4vn_vce = NULL) {
  sig <- .r4vn_lr_signature(model, label)

  ll_obj <- tryCatch(stats::logLik(model), error = function(e) NULL)
  if (is.null(ll_obj)) {
    stop(
      sprintf("A log likelihood could not be obtained from `%s`.", label),
      call. = FALSE
    )
  }

  ll <- as.numeric(ll_obj)[1L]
  if (!is.finite(ll)) {
    stop(
      sprintf("`%s` does not have a finite log likelihood.", label),
      call. = FALSE
    )
  }

  parameters <- attr(ll_obj, "df")
  if (is.null(parameters) || !length(parameters) || !is.finite(parameters)) {
    parameters <- sum(!is.na(stats::coef(model)))
  }
  parameters <- as.numeric(parameters)[1L]

  f <- tryCatch(stats::formula(model), error = function(e) NULL)
  ftext <- if (is.null(f)) "" else .r4vn_deparse1(f)

  term_labels <- tryCatch(
    attr(stats::terms(model), "term.labels"),
    error = function(e) character(0)
  )

  list(
    model = model,
    label = label,
    type = sig$type,
    signature = sig$signature,
    logLik = ll,
    parameters = parameters,
    n = as.numeric(stats::nobs(model)),
    AIC = tryCatch(as.numeric(stats::AIC(model))[1L], error = function(e) NA_real_),
    BIC = tryCatch(as.numeric(stats::BIC(model))[1L], error = function(e) NA_real_),
    formula = ftext,
    terms = as.character(term_labels),
    model.frame = .r4vn_lr_model_frame(model),
    model.matrix = .r4vn_lr_model_matrix(model),
    vce = r4vn_vce
  )
}

.r4vn_lr_same_vector <- function(x, y) {
  isTRUE(all.equal(
    x, y,
    check.attributes = TRUE
  ))
}

.r4vn_lr_is_nested <- function(reduced, full) {
  xr <- reduced$model.matrix
  xf <- full$model.matrix

  if (!is.null(xr) && !is.null(xf) &&
      nrow(xr) == nrow(xf) &&
      ncol(xr) > 0L && ncol(xf) > 0L &&
      all(is.finite(xr)) && all(is.finite(xf))) {
    rank_full <- qr(xf, tol = 1e-9)$rank
    rank_aug <- qr(cbind(xf, xr), tol = 1e-9)$rank
    return(identical(rank_full, rank_aug))
  }

  all(reduced$terms %in% full$terms)
}

.r4vn_lr_compare_pair <- function(reduced, full,
                                  reduced_label, full_label) {
  if (!identical(reduced$signature, full$signature)) {
    stop(
      sprintf(
        "`%s` and `%s` are not the same model type/family/link.",
        reduced_label, full_label
      ),
      call. = FALSE
    )
  }

  if (!identical(as.numeric(reduced$n), as.numeric(full$n))) {
    stop(
      sprintf(
        paste0(
          "`%s` and `%s` were fitted to different observations ",
          "(n = %s versus n = %s). Refit both models on the same analytic sample."
        ),
        reduced_label, full_label, reduced$n, full$n
      ),
      call. = FALSE
    )
  }

  mf1 <- reduced$model.frame
  mf2 <- full$model.frame

  if (!is.null(mf1) && !is.null(mf2)) {
    if (!identical(row.names(mf1), row.names(mf2))) {
      stop(
        sprintf(
          "`%s` and `%s` contain different observations. Refit both models on the same analytic sample.",
          reduced_label, full_label
        ),
        call. = FALSE
      )
    }

    y1 <- tryCatch(stats::model.response(mf1), error = function(e) NULL)
    y2 <- tryCatch(stats::model.response(mf2), error = function(e) NULL)
    if (!is.null(y1) && !is.null(y2) && !.r4vn_lr_same_vector(y1, y2)) {
      stop(
        sprintf(
          "`%s` and `%s` do not use the same outcome.",
          reduced_label, full_label
        ),
        call. = FALSE
      )
    }

    w1 <- stats::model.weights(mf1)
    w2 <- stats::model.weights(mf2)
    if (is.null(w1)) w1 <- rep(1, nrow(mf1))
    if (is.null(w2)) w2 <- rep(1, nrow(mf2))
    if (!.r4vn_lr_same_vector(as.numeric(w1), as.numeric(w2))) {
      stop(
        sprintf(
          "`%s` and `%s` use different model weights.",
          reduced_label, full_label
        ),
        call. = FALSE
      )
    }

    o1 <- stats::model.offset(mf1)
    o2 <- stats::model.offset(mf2)
    if (is.null(o1)) o1 <- rep(0, nrow(mf1))
    if (is.null(o2)) o2 <- rep(0, nrow(mf2))
    if (!.r4vn_lr_same_vector(as.numeric(o1), as.numeric(o2))) {
      stop(
        sprintf(
          "`%s` and `%s` use different offsets/exposure definitions.",
          reduced_label, full_label
        ),
        call. = FALSE
      )
    }
  }

  if (full$parameters <= reduced$parameters) {
    stop(
      sprintf(
        paste0(
          "`%s` must be the reduced model and `%s` the larger nested model. ",
          "Supply models from smallest to largest."
        ),
        reduced_label, full_label
      ),
      call. = FALSE
    )
  }

  if (!.r4vn_lr_is_nested(reduced, full)) {
    stop(
      sprintf(
        "`%s` is not nested within `%s`.",
        reduced_label, full_label
      ),
      call. = FALSE
    )
  }

  statistic <- 2 * (full$logLik - reduced$logLik)

  if (statistic < -1e-7) {
    stop(
      sprintf(
        paste0(
          "The larger model `%s` has a lower log likelihood than `%s`. ",
          "Check model convergence and nesting."
        ),
        full_label, reduced_label
      ),
      call. = FALSE
    )
  }

  statistic <- max(statistic, 0)
  df <- as.integer(round(full$parameters - reduced$parameters))

  if (df < 1L) {
    stop(
      "The difference in model degrees of freedom must be positive.",
      call. = FALSE
    )
  }

  p <- stats::pchisq(
    statistic,
    df = df,
    lower.tail = FALSE
  )

  new_terms <- setdiff(full$terms, reduced$terms)
  added <- if (length(new_terms)) {
    paste(.r4vn_mx_clean_term(new_terms), collapse = ", ")
  } else {
    "Additional parameters"
  }

  list(
    statistic = statistic,
    df = df,
    p.value = p,
    added = added
  )
}


# ============================================================================
# Function source: lrtest.R (consolidated here)
# ============================================================================
#' Likelihood-ratio test for nested regression models
#'
#' Compares two or more nested likelihood-based regression models. Objects
#' returned by R4VN [logistic()] and [poisson()] can be supplied directly.
#'
#' @param ... Two or more nested fitted models, ordered from the smaller model
#'   to progressively larger models. R4VN statistical results containing a
#'   fitted model in `raw$model` are accepted directly.
#' @param digits Number of decimal places for likelihood and LR statistics.
#' @param p_digits Number of decimal places for p-values.
#' @param show Logical; open the formatted result in the Viewer. Default `TRUE`.
#' @param console Logical; also print the result in the Console. Default `FALSE`.
#'
#' @details
#' When more than two models are supplied, comparisons are sequential:
#' `M1` versus `M2`, then `M2` versus `M3`, and so on.
#'
#' The models must use the same outcome, analytic observations, weights,
#' offsets/exposure definition, and likelihood family/link, and each larger
#' model must contain the smaller model.
#'
#' Supported fits include ordinary likelihood-based `glm` models such as
#' logistic and Poisson regression, `MASS::glm.nb()` negative-binomial models,
#' `survival::coxph()` Cox models, and `survival::survreg()` parametric survival
#' models.
#'
#' Quasi-likelihood models are not supported. R4VN models fitted with
#' `vce = "robust"` or `vce = "cluster"` are also rejected because the
#' classical likelihood-ratio chi-square test is model-likelihood inference,
#' not robust covariance inference.
#'
#' For ordinary linear regression use the nested-model F test rather than
#' `lrtest()`.
#'
#' @return An object of class `r4vn_stat`. The unformatted comparison table is
#'   stored in `result$raw$table`; the backward-compatible `result$raw$comparison` table is also retained.
#'
#' @examples
#' set.seed(2026)
#' d <- data.frame(
#'   y = factor(rbinom(250, 1, .35), levels = 0:1,
#'              labels = c("No", "Yes")),
#'   age = rnorm(250, 45, 12),
#'   sex = factor(sample(c("Female", "Male"), 250, TRUE)),
#'   treatment = factor(sample(c("No", "Yes"), 250, TRUE))
#' )
#'
#' m1 <- logistic(y, c.age, i.sex, i.treatment,
#'                data = d, event = "Yes", show = FALSE)
#' m2 <- logistic(y, c.age, i.sex*i.treatment,
#'                data = d, event = "Yes", show = FALSE)
#'
#' lrtest(m1, m2, show = FALSE)
#'
#' @seealso [logistic()], [poisson()]
#' @export
lrtest <- function(...,
                   digits = 3,
                   p_digits = 3,
                   show = TRUE,
                   console = FALSE) {
  exprs <- as.list(substitute(list(...)))[-1L]
  objects <- list(...)

  if (length(objects) < 2L) {
    stop("`lrtest()` requires at least two models.", call. = FALSE)
  }

  expr_names <- names(exprs)
  if (is.null(expr_names)) {
    expr_names <- rep("", length(exprs))
  }

  labels <- vapply(
    seq_along(exprs),
    function(i) {
      if (nzchar(expr_names[i])) {
        expr_names[i]
      } else {
        .r4vn_deparse1(exprs[[i]])
      }
    },
    character(1)
  )

  extracted <- lapply(
    seq_along(objects),
    function(i) .r4vn_lr_extract(objects[[i]], labels[i])
  )

  info <- lapply(
    seq_along(extracted),
    function(i) {
      .r4vn_lr_model_info(
        extracted[[i]]$model,
        labels[i],
        extracted[[i]]$vce
      )
    }
  )

  signatures <- vapply(
    info,
    function(z) z$signature,
    character(1)
  )

  if (length(unique(signatures)) != 1L) {
    stop(
      "The supplied models are not of the same likelihood-based model type, distribution, and link.",
      call. = FALSE
    )
  }

  k <- length(info)
  lr_stat <- rep(NA_real_, k)
  lr_df <- rep(NA_integer_, k)
  lr_p <- rep(NA_real_, k)
  added <- rep("", k)

  for (i in 2:k) {
    cmp <- .r4vn_lr_compare_pair(
      info[[i - 1L]],
      info[[i]],
      labels[i - 1L],
      labels[i]
    )

    lr_stat[i] <- cmp$statistic
    lr_df[i] <- cmp$df
    lr_p[i] <- cmp$p.value
    added[i] <- cmp$added
  }

  # Current consolidated/raw contract.
  raw_table <- data.frame(
    Model = paste0("M", seq_len(k)),
    Object = labels,
    n = vapply(info, function(z) z$n, numeric(1)),
    Parameters = vapply(info, function(z) z$parameters, numeric(1)),
    logLik = vapply(info, function(z) z$logLik, numeric(1)),
    AIC = vapply(info, function(z) z$AIC, numeric(1)),
    BIC = vapply(info, function(z) z$BIC, numeric(1)),
    Added = added,
    LR = lr_stat,
    df = lr_df,
    p = lr_p,
    stringsAsFactors = FALSE,
    check.names = FALSE
  )

  # Backward-compatible raw contract used by the earlier R4VN LR update and
  # Studio integrations. Both are intentionally retained.
  raw_comparison <- data.frame(
    Model = raw_table$Model,
    Terms.added = ifelse(seq_len(k) == 1L, "-", added),
    Log.likelihood = raw_table$logLik,
    df = raw_table$Parameters,
    LR.chi2 = raw_table$LR,
    LR.df = raw_table$df,
    p.value = raw_table$p,
    AIC = raw_table$AIC,
    BIC = raw_table$BIC,
    stringsAsFactors = FALSE,
    check.names = FALSE
  )

  lr_show <- rep("", k)
  df_show <- rep("", k)
  p_show <- rep("", k)

  if (k >= 2L) {
    lr_show[2:k] <- .r4vn_num(lr_stat[2:k], digits)
    df_show[2:k] <- as.character(lr_df[2:k])
    p_show[2:k] <- .r4vn_p(lr_p[2:k], p_digits)
  }

  added_show <- added
  added_show[1L] <- "-"

  display_table <- data.frame(
    Model = paste0("M", seq_len(k)),
    `Added terms` = added_show,
    n = raw_table$n,
    `Log likelihood` = .r4vn_num(raw_table$logLik, digits),
    AIC = .r4vn_num(raw_table$AIC, digits),
    `LR chi2` = lr_show,
    df = df_show,
    p = p_show,
    stringsAsFactors = FALSE,
    check.names = FALSE
  )

  formula_table <- data.frame(
    Model = paste0("M", seq_len(k)),
    Object = labels,
    Formula = vapply(
      info,
      function(z) z$formula,
      character(1)
    ),
    stringsAsFactors = FALSE,
    check.names = FALSE
  )

  note <- paste0(
    "Each model after M1 is compared with the immediately preceding model. ",
    "Models must be nested and fitted to the same analytic observations."
  )

  if (identical(info[[1L]]$type, "Cox regression")) {
    note <- paste0(
      note,
      " Cox comparisons use the partial likelihood."
    )
  }

  result <- .r4vn_result(
    "Likelihood-ratio test",
    sections = list(
      "Model comparison" = display_table,
      "Models" = formula_table
    ),
    notes = note,
    raw = list(
      table = raw_table,
      comparison = raw_comparison,
      models = lapply(extracted, function(z) z$model),
      type = info[[1L]]$type
    ),
    call = match.call()
  )

  .r4vn_show(
    result,
    show = show,
    console = console
  )
}


# ============================================================================
# Viewer renderers for model commands
# ============================================================================
.r4vn_model_viewer <- function(x, subtitle) {
  blocks <- character()
  nms <- names(x$sections)
  if ("Model summary" %in% nms) {
    blocks <- c(blocks, paste0('<section class="r4vn-section"><h2>Model summary</h2>',
                               .r4vn_view_key_values(x$sections[["Model summary"]]), '</section>'))
    nms <- setdiff(nms, "Model summary")
  }
  for (nm in nms) blocks <- c(blocks, .r4vn_view_section(nm, x$sections[[nm]]))
  .r4vn_view_document(x$title, paste0(blocks, collapse = ""), notes = x$notes,
                      subtitle = subtitle, prefix = "r4vn-model-")
}

.r4vn_viewer_corr <- function(x) {
  blocks <- character()
  if (!is.null(x$sections[["Correlation matrix"]]))
    blocks <- c(blocks, .r4vn_view_section("Correlation matrix", x$sections[["Correlation matrix"]]))
  if (!is.null(x$sections[["P-values"]]))
    blocks <- c(blocks, .r4vn_view_section("P-values", x$sections[["P-values"]]))
  if (!is.null(x$sections[["Observations"]]))
    blocks <- c(blocks, .r4vn_view_section("Pairwise observations", x$sections[["Observations"]]))
  if (!is.null(x$sections[["Pairwise confidence intervals"]]))
    blocks <- c(blocks, .r4vn_view_section("Pairwise confidence intervals", x$sections[["Pairwise confidence intervals"]]))
  .r4vn_view_document(x$title, paste0(blocks, collapse = ""), notes = x$notes,
                      subtitle = "Correlation analysis", prefix = "r4vn-corr-")
}

.r4vn_viewer_regress <- function(x) .r4vn_model_viewer(x, "Linear regression model")
.r4vn_viewer_logistic <- function(x) .r4vn_model_viewer(x, "Binary logistic regression model")
.r4vn_viewer_poisson <- function(x) .r4vn_model_viewer(x, "Poisson regression model")
.r4vn_viewer_lrtest <- function(x) .r4vn_model_viewer(x, "Nested likelihood-based model comparison")