R/postestimation.R

In R4VN: Health Data Analysis and Publication-Ready Reporting

# ==========================================================================
# Postestimation commands: at(), margins(), marginsplot(), predict(), lincom()
# ==========================================================================

#' Specify predictor values for margins
#' @usage at(...)
#'
#' @description
#' `at()` is a small helper used by [margins()]. Each argument names a predictor
#' and supplies one or more values. Several vectors are crossed automatically.
#'
#' @param ... Named predictor values, for example `at(age = 40)` or
#'   `at(age = seq(30, 50, 5), hypertension = c(0, 1))`.
#' @return An object of class `r4vn_at`.
#' @export
#' @examples
#' at(age = 40)
#' at(age = seq(30, 50, 5), hypertension = c(0, 1))
at <- function(...) {
  z <- list(...)
  if (!length(z) || is.null(names(z)) || any(!nzchar(names(z)))) {
    stop("Every `at()` value must be named, for example `at(age = 40)`.", call. = FALSE)
  }
  if (anyDuplicated(names(z))) stop("Predictors in `at()` must be unique.", call. = FALSE)
  if (any(lengths(z) == 0L)) stop("Every `at()` predictor must have at least one value.", call. = FALSE)
  structure(z, class = c("r4vn_at", "list"))
}

.r4vn_unwrap_model <- function(model = NULL, env = parent.frame(), expr = substitute(model)) {
  if (missing(model) || is.null(model)) return(.r4vn_get_active_model(TRUE))
  obj <- model
  if (inherits(obj, "r4vn_stat") && !is.null(obj$raw$model)) obj <- obj$raw$model
  if (is.list(obj) && !inherits(obj, c("lm", "glm", "coxph", "rq", "aov", "polr", "multinom")) && !is.null(obj$model) && inherits(obj$model, c("lm", "glm", "coxph", "rq", "aov"))) obj <- obj$model
  if (is.null(obj)) stop("A fitted model is required.", call. = FALSE)
  obj
}

.r4vn_post_vcov <- function(model) {
  active <- .r4vn_get_active_model(required = FALSE)
  res <- .r4vn_get_active_model(required = FALSE, result = TRUE)
  if (!is.null(active) && identical(model, active) && is.list(res) && is.list(res$raw) && !is.null(res$raw$vcov)) {
    return(as.matrix(res$raw$vcov))
  }
  tryCatch(as.matrix(stats::vcov(model)), error = function(e) NULL)
}

.r4vn_model_base_data <- function(model) {
  mf <- tryCatch(stats::model.frame(model), error = function(e) NULL)
  tt <- tryCatch(stats::terms(model), error = function(e) NULL)
  raw_pred <- tryCatch(model$.r4vn_prediction_data, error = function(e) NULL)
  if (is.data.frame(raw_pred) && nrow(raw_pred)) {
    response <- if (!is.null(tt) && attr(tt, "response") > 0L && !is.null(mf)) names(mf)[attr(tt, "response")] else NULL
    return(list(data = raw_pred, response = response, frame = mf, terms = tt))
  }
  if (is.null(mf) || !is.data.frame(mf)) stop("The fitted model does not retain a usable model frame.", call. = FALSE)
  response <- if (!is.null(tt) && attr(tt, "response") > 0L) names(mf)[attr(tt, "response")] else NULL
  pred <- mf
  if (!is.null(response) && response %in% names(pred)) pred[[response]] <- NULL
  list(data = pred, response = response, frame = mf, terms = tt)
}

# Coerce user/new active data to the predictor storage types retained by the
# fitted model. R4VN compact syntax can deliberately turn a numeric 0/1 column
# into a factor (for example i.htn) without mutating the user's data frame.
# Base predict.glm() correctly rejects numeric newdata for a factor term, but
# that is surprising in a Stata-like postestimation workflow. Reconstruct the
# fitted predictor types here while preserving row count/order and rejecting
# genuinely unknown levels.
.r4vn_prepare_prediction_data <- function(model, newdata) {
  if (!is.data.frame(newdata)) return(newdata)
  ref <- tryCatch(.r4vn_model_base_data(model)$data, error = function(e) NULL)
  if (is.null(ref) || !is.data.frame(ref)) return(newdata)

  out <- newdata

  # The fitted model's xlevels are the authoritative record of predictors that
  # were categorical at estimation time. This is essential for R4VN compact
  # syntax such as i.htn: the user's source data can remain numeric 0/1 while
  # the fitted model legitimately expects a factor. Rebuild those factor
  # columns before calling the underlying predict method.
  xlev <- tryCatch(model$xlevels, error = function(e) NULL)
  factor_predictors <- character()
  if (is.list(xlev) && length(xlev)) {
    for (term_name in names(xlev)) {
      nm <- term_name
      # Usually xlevels is named with the original variable (e.g. htn). This
      # fallback also handles ordinary formulas written as factor(htn).
      if (!nm %in% names(out) && grepl("^factor\\([^()]+\\)$", nm)) {
        nm2 <- sub("^factor\\(([^()]+)\\)$", "\\1", nm)
        if (nm2 %in% names(out)) nm <- nm2
      }
      if (!nm %in% names(out)) next
      allowed <- as.character(xlev[[term_name]])
      vals <- as.character(out[[nm]])
      bad <- unique(vals[!is.na(vals) & !vals %in% allowed])
      if (length(bad)) {
        stop(
          "New data contain unknown level(s) for factor `", nm, "`: ",
          paste(bad, collapse = ", "), ". Allowed levels: ",
          paste(allowed, collapse = ", "), ".",
          call. = FALSE
        )
      }
      out[[nm]] <- factor(vals, levels = allowed)
      factor_predictors <- unique(c(factor_predictors, nm))
    }
  }

  for (nm in intersect(names(ref), names(out))) {
    template <- ref[[nm]]
    x <- out[[nm]]

    # Do not undo the authoritative factor conversion above merely because a
    # secondary prediction-data cache retained the original numeric storage.
    if (nm %in% factor_predictors) next

    if (is.factor(template)) {
      vals <- as.character(x)
      allowed <- levels(template)
      bad <- unique(vals[!is.na(vals) & !vals %in% allowed])
      if (length(bad)) {
        stop(
          "New data contain unknown level(s) for factor `", nm, "`: ",
          paste(bad, collapse = ", "), ". Allowed levels: ",
          paste(allowed, collapse = ", "), ".",
          call. = FALSE
        )
      }
      out[[nm]] <- factor(vals, levels = allowed, ordered = is.ordered(template))
      next
    }

    if (is.logical(template) && !is.logical(x)) {
      vals <- if (is.numeric(x)) {
        ifelse(is.na(x), NA, x != 0)
      } else {
        xx <- tolower(trimws(as.character(x)))
        ans <- rep(NA, length(xx))
        ans[xx %in% c("true", "t", "1", "yes", "y")] <- TRUE
        ans[xx %in% c("false", "f", "0", "no", "n")] <- FALSE
        ans
      }
      out[[nm]] <- as.logical(vals)
      next
    }

    if (is.numeric(template) && !is.numeric(x)) {
      vals <- suppressWarnings(as.numeric(as.character(x)))
      bad <- !is.na(x) & is.na(vals)
      if (any(bad)) {
        stop("New data for numeric predictor `", nm, "` could not be converted to numeric.", call. = FALSE)
      }
      out[[nm]] <- vals
      next
    }

    # Preserve Date/POSIX classes where the model frame retained them.
    if (inherits(template, "Date") && !inherits(x, "Date")) {
      out[[nm]] <- as.Date(x)
    } else if (inherits(template, "POSIXct") && !inherits(x, "POSIXct")) {
      out[[nm]] <- as.POSIXct(x)
    }
  }

  # Cox models created by tabsurv()/cox() use safe internal column names. Add
  # those columns after coercing the original user-facing predictors above.
  pmap <- tryCatch(model$.r4vn_prediction_map, error = function(e) NULL)
  if (is.data.frame(pmap) && nrow(pmap) && all(c("variable", "internal") %in% names(pmap))) {
    mf <- tryCatch(stats::model.frame(model), error = function(e) NULL)
    for (i in seq_len(nrow(pmap))) {
      src <- as.character(pmap$variable[i])
      dst <- as.character(pmap$internal[i])
      if (!src %in% names(out)) next
      z <- out[[src]]
      template <- if (is.data.frame(mf) && dst %in% names(mf)) mf[[dst]] else NULL
      if (is.factor(template)) {
        vals <- as.character(z)
        allowed <- levels(template)
        bad <- unique(vals[!is.na(vals) & !vals %in% allowed])
        if (length(bad)) {
          stop("New data contain unknown level(s) for factor `", src, "`: ",
               paste(bad, collapse = ", "), ". Allowed levels: ",
               paste(allowed, collapse = ", "), ".", call. = FALSE)
        }
        z <- factor(vals, levels = allowed, ordered = is.ordered(template))
      } else if (is.numeric(template) && !is.numeric(z)) {
        zz <- suppressWarnings(as.numeric(as.character(z)))
        if (any(!is.na(z) & is.na(zz))) {
          stop("New data for numeric predictor `", src, "` could not be converted to numeric.", call. = FALSE)
        }
        z <- zz
      }
      out[[dst]] <- z
    }
  }
  out
}

.r4vn_at_grid <- function(at_spec, model_data) {
  if (is.null(at_spec)) return(data.frame(.r4vn_scenario = 1L)[FALSE, , drop = FALSE])
  if (!inherits(at_spec, "r4vn_at") && !is.list(at_spec)) stop("`at` must be created by `at()` or be a named list.", call. = FALSE)
  if (is.null(names(at_spec)) || any(!nzchar(names(at_spec)))) stop("Every `at` predictor must be named.", call. = FALSE)
  miss <- setdiff(names(at_spec), names(model_data))
  if (length(miss)) stop("Predictor(s) in `at` were not found in the fitted model: ", paste(miss, collapse = ", "), call. = FALSE)
  vals <- lapply(names(at_spec), function(nm) {
    v <- at_spec[[nm]]
    base <- model_data[[nm]]
    if (is.factor(base)) {
      vv <- as.character(v)
      bad <- setdiff(vv, levels(base))
      if (length(bad)) stop("Unknown factor level(s) for `", nm, "`: ", paste(bad, collapse = ", "), call. = FALSE)
      factor(vv, levels = levels(base), ordered = is.ordered(base))
    } else v
  })
  names(vals) <- names(at_spec)
  do.call(expand.grid, c(vals, list(KEEP.OUT.ATTRS = FALSE, stringsAsFactors = FALSE)))
}

.r4vn_mode_value <- function(x) {
  z <- x[!is.na(x)]
  if (!length(z)) return(NA)
  tb <- table(z, useNA = "no"); names(tb)[which.max(tb)]
}

.r4vn_representative_row <- function(d) {
  out <- d[1L, , drop = FALSE]
  for (nm in names(d)) {
    z <- d[[nm]]
    if (is.numeric(z)) out[[nm]] <- mean(z, na.rm = TRUE)
    else if (is.factor(z)) out[[nm]] <- factor(.r4vn_mode_value(z), levels = levels(z), ordered = is.ordered(z))
    else if (is.logical(z)) out[[nm]] <- as.logical(.r4vn_mode_value(z))
    else out[[nm]] <- .r4vn_mode_value(z)
  }
  out
}

.r4vn_design_matrix <- function(model, newdata) {
  tt <- tryCatch(stats::delete.response(stats::terms(model)), error = function(e) NULL)
  if (is.null(tt)) return(NULL)
  newdata <- .r4vn_prepare_prediction_data(model, newdata)
  tryCatch(stats::model.matrix(tt, newdata, contrasts.arg = model$contrasts, xlev = model$xlevels), error = function(e) NULL)
}

.r4vn_margin_one <- function(model, nd, type = c("response", "link"), level = 0.95, V = NULL) {
  type <- match.arg(type)
  nd_model <- .r4vn_prepare_prediction_data(model, nd)
  is_glm <- inherits(model, "glm")
  is_cox <- inherits(model, "coxph")
  is_lm <- inherits(model, "lm") && !is_glm

  if (is_cox) {
    lp <- tryCatch(as.numeric(stats::predict(model, newdata = nd_model, type = "lp", reference = "zero")), error = function(e) NULL)
    if (is.null(lp)) stop("Predictions could not be calculated for a margins scenario.", call. = FALSE)
    pred <- if (type == "response") exp(lp) else lp
  } else if (is_glm) {
    pred <- tryCatch(as.numeric(stats::predict(model, newdata = nd_model, type = type)), error = function(e) NULL)
    if (is.null(pred)) stop("Predictions could not be calculated for a margins scenario.", call. = FALSE)
  } else {
    # lm/aov/rq and other scalar-prediction models generally use their natural
    # response scale and may not implement a GLM-like `type` argument.
    pred <- tryCatch(as.numeric(stats::predict(model, newdata = nd_model)), error = function(e) NULL)
    if (is.null(pred)) stop("Predictions could not be calculated for a margins scenario.", call. = FALSE)
  }

  est <- mean(pred, na.rm = TRUE)
  se <- NA_real_
  X <- .r4vn_design_matrix(model, nd_model)
  beta <- tryCatch(stats::coef(model), error = function(e) NULL)
  if (!is.null(X) && !is.null(beta) && !is.null(V)) {
    common <- intersect(intersect(colnames(X), names(beta)), intersect(rownames(V), colnames(V)))
    common <- common[is.finite(beta[common])]
    if (length(common)) {
      XX <- X[, common, drop = FALSE]
      VV <- V[common, common, drop = FALSE]
      if (is_glm && type == "response") {
        eta <- tryCatch(as.numeric(stats::predict(model, newdata = nd_model, type = "link")), error = function(e) NULL)
        grad <- if (is.null(eta)) NULL else colMeans(XX * stats::family(model)$mu.eta(eta), na.rm = TRUE)
      } else if (is_cox && type == "response") {
        # d exp(Xb) / db = exp(Xb) X. reference="zero" above prevents
        # prediction centering from contaminating this derivative.
        grad <- colMeans(XX * pred, na.rm = TRUE)
      } else {
        grad <- colMeans(XX, na.rm = TRUE)
      }
      if (!is.null(grad) && all(is.finite(grad))) {
        vv <- drop(t(grad) %*% VV %*% grad)
        if (is.finite(vv) && vv >= 0) se <- sqrt(vv)
      }
    }
  }

  crit <- if (is_lm && is.finite(stats::df.residual(model)) && stats::df.residual(model) > 0) {
    stats::qt(1 - (1 - level) / 2, df = stats::df.residual(model))
  } else stats::qnorm(1 - (1 - level) / 2)
  ci <- est + c(-1, 1) * crit * se
  if (is_glm && type == "response") {
    fam <- stats::family(model)$family
    if (identical(fam, "binomial")) ci <- pmax(0, pmin(1, ci))
    if (identical(fam, "poisson")) ci[1L] <- max(0, ci[1L])
  }
  if (is_cox && type == "response") ci[1L] <- max(0, ci[1L])
  c(estimate = est, se = se, lower = ci[1L], upper = ci[2L], n = sum(is.finite(pred)))
}

.r4vn_resolve_over <- function(expr, model_data, env) {
  if (.r4vn_expr_is_null(expr)) return(character())
  .r4vn_resolve_name_spec(expr, model_data, env, "over", allow_null = TRUE, multiple = TRUE)
}

#' Predictive margins after an R4VN model
#' @usage
#' margins(
#'   model = NULL, at = NULL, over = NULL, type = c("response", "link"),
#'   atmeans = FALSE, level = 0.95, digits = 3, p_digits = 3, show = TRUE,
#'   console = FALSE
#' )
#'
#' @description
#' Calculates average adjusted predictions from the active (most recently fitted)
#' model or from a model supplied explicitly. This follows the predictive-margin
#' interpretation used by Stata: covariates not listed in `at()` retain their
#' observed values and predictions are averaged over the estimation sample.
#'
#' @param model Optional fitted model or R4VN model result. When omitted, the
#'   most recently fitted R4VN model is used.
#' @param at Predictor settings created with [at()]. Vectors are crossed to form
#'   all requested scenarios.
#' @param over Optional grouping variable(s). `vars(a, b)` produces margins for
#'   observed combinations of `a` and `b`.
#' @param type `"response"` for adjusted outcomes/probabilities/rates (relative risk for Cox models) or `"link"`
#'   for the linear predictor.
#' @param atmeans If `TRUE`, covariates not fixed in `at()` are set to their mean
#'   (numeric) or modal category rather than averaged over individual observations.
#' @param level Confidence level.
#' @param digits,p_digits Formatting controls.
#' @param show,console Display controls.
#' @return An object of class `r4vn_margins` and `r4vn_stat`.
#' @export
#' @examples
#' d <- data.frame(
#'   outcome = factor(c(0, 0, 0, 1, 0, 1, 1, 1), levels = 0:1,
#'                    labels = c("No", "Yes")),
#'   age = c(20, 25, 30, 35, 40, 45, 50, 55),
#'   sex = factor(rep(c("Female", "Male"), 4))
#' )
#' m <- logistic(outcome, c.age, i.sex, data = d, event = "Yes", show = FALSE)
#' margins(m, at = at(age = 40), show = FALSE)
#' margins(m, over = sex, show = FALSE)
margins <- function(model = NULL, at = NULL, over = NULL, type = c("response", "link"),
                    atmeans = FALSE, level = 0.95, digits = 3, p_digits = 3,
                    show = TRUE, console = FALSE) {
  call <- match.call(); env <- parent.frame(); type <- match.arg(type)
  if (!is.numeric(level) || length(level) != 1L || level <= 0 || level >= 1) stop("`level` must be between 0 and 1.", call. = FALSE)
  fit <- if (missing(model) || is.null(model)) .r4vn_get_active_model(TRUE) else .r4vn_unwrap_model(model)
  base <- .r4vn_model_base_data(fit); md <- base$data
  at_spec <- if (missing(at)) NULL else at
  grid <- .r4vn_at_grid(at_spec, md)
  if (!nrow(grid)) grid <- data.frame(.r4vn_one = 1L)[, FALSE, drop = FALSE]
  over_names <- .r4vn_resolve_over(substitute(over), md, env)
  at_names0 <- if (is.null(at_spec)) character() else names(at_spec)
  overlap <- intersect(at_names0, over_names)
  if (length(overlap)) {
    stop("A predictor cannot be specified in both `at` and `over`: ", paste(overlap, collapse = ", "), ".", call. = FALSE)
  }
  over_grid <- if (!length(over_names)) data.frame(.r4vn_one = 1L)[, FALSE, drop = FALSE] else unique(md[stats::complete.cases(md[over_names]), over_names, drop = FALSE])
  if (!nrow(over_grid)) over_grid <- data.frame(.r4vn_one = 1L)[, FALSE, drop = FALSE]
  scen <- merge(grid, over_grid, by = NULL)
  if (!nrow(scen)) scen <- data.frame(.r4vn_one = 1L)[, FALSE, drop = FALSE]
  V <- .r4vn_post_vcov(fit)
  rows <- vector("list", nrow(scen))

  for (i in seq_len(nrow(scen))) {
    nd <- md
    if (length(over_names)) {
      keep <- rep(TRUE, nrow(nd))
      for (nm in over_names) keep <- keep & !is.na(nd[[nm]]) & as.character(nd[[nm]]) == as.character(scen[[nm]][i])
      nd <- nd[keep, , drop = FALSE]
    }
    if (!nrow(nd)) next
    if (isTRUE(atmeans)) nd <- .r4vn_representative_row(nd)
    at_names <- if (is.null(at_spec)) character() else names(at_spec)
    for (nm in at_names) {
      value <- scen[[nm]][i]
      if (is.factor(md[[nm]])) nd[[nm]] <- factor(rep(as.character(value), nrow(nd)), levels = levels(md[[nm]]), ordered = is.ordered(md[[nm]]))
      else nd[[nm]] <- rep(value, nrow(nd))
    }
    z <- .r4vn_margin_one(fit, nd, type = type, level = level, V = V)
    prefix <- if (ncol(scen)) scen[i, , drop = FALSE] else data.frame()
    rows[[i]] <- cbind(prefix, data.frame(Margin = unname(z["estimate"]), SE = unname(z["se"]), Lower = unname(z["lower"]), Upper = unname(z["upper"]), N = unname(z["n"]), stringsAsFactors = FALSE))
  }
  rows <- Filter(Negate(is.null), rows)
  if (!length(rows)) stop("No margins were estimable for the requested scenarios.", call. = FALSE)
  tab_raw <- do.call(rbind, rows); rownames(tab_raw) <- NULL
  tab <- tab_raw
  for (nm in intersect(c("Margin", "SE", "Lower", "Upper"), names(tab))) tab[[nm]] <- .r4vn_num(tab[[nm]], digits)
  notes <- c(
    "Margins are average predictions over the estimation sample unless atmeans = TRUE.",
    if (length(over_names)) paste0("Margins are reported within observed `over` groups: ", paste(over_names, collapse = ", "), ".") else NULL
  )
  out <- .r4vn_result("Predictive margins", list("Margins" = tab), notes = notes,
    raw = list(model = fit, margins = tab_raw, at = at_spec, over = over_names, type = type, atmeans = atmeans, vcov = V), call = call)
  class(out) <- unique(c("r4vn_margins", class(out)))
  .r4vn_set_last_margins(out)
  .r4vn_show(out, show = show, console = console)
}

#' Plot predictive margins
#' @usage
#' marginsplot(
#'   result = NULL, x = NULL, by = NULL, ci = TRUE, line = TRUE, points = TRUE,
#'   line_width = 2, line_type = 1, point_size = 1, point_pch = 16,
#'   ci_color = NULL, ci_alpha = 1, ci_lwd = 1, ci_lty = 1, xline = NULL,
#'   yline = NULL, ref_color = "gray40", ref_lty = 2, ref_lwd = 1, xlab = NULL,
#'   ylab = NULL, xtitle = NULL, ytitle = NULL, title = NULL, subtitle = NULL,
#'   note = NULL, color = NULL, palette = "journal", alpha = 1, legend = TRUE,
#'   legend_position = "topright", theme = "journal", size = 11, file = NULL,
#'   width = 7, height = 5, dpi = 300, show = TRUE, bg = "white", vline = NULL,
#'   hline = NULL
#' )
#'
#' @param result A result returned by [margins()]. When omitted, the most recent
#'   margins result is used.
#' @param x Optional x-axis variable from the margins table. The first varying
#'   `at()` variable is chosen automatically.
#' @param by Optional grouping variable from the margins table.
#' @param ci Draw confidence-interval error bars.
#' @param line,points Draw connecting lines and points.
#' @param line_width,line_type Connecting-line width and line type.
#' @param point_size,point_pch Point size and symbol.
#' @param ci_color,ci_alpha,ci_lwd,ci_lty Confidence-interval colour, transparency, width, and line type. `ci_color = NULL` follows each group colour.
#' @param xline,yline Optional reference lines at x- and y-axis values.
#' @param vline,hline Deprecated aliases for `xline` and `yline`.
#' @param ref_color,ref_lty,ref_lwd Reference-line formatting.
#' @param xlab,ylab,xtitle,ytitle,title,subtitle,note,color,palette,alpha,legend,legend_position,theme,size,file,width,height,dpi,show,bg Graph controls consistent with other R4VN graphs.
#' @return An `r4vn_graph` object.
#' @examples
#' d <- data.frame(
#'   outcome = factor(c(0, 0, 0, 1, 0, 1, 1, 1), levels = 0:1,
#'                    labels = c("No", "Yes")),
#'   age = c(20, 25, 30, 35, 40, 45, 50, 55)
#' )
#' m <- logistic(outcome, c.age, data = d, event = "Yes", show = FALSE)
#' mg <- margins(m, at = at(age = seq(30, 50, 10)), show = FALSE)
#' marginsplot(mg, show = FALSE)
#' @export
marginsplot <- function(result = NULL, x = NULL, by = NULL, ci = TRUE, line = TRUE, points = TRUE,
                        line_width = 2, line_type = 1, point_size = 1, point_pch = 16,
                        ci_color = NULL, ci_alpha = 1, ci_lwd = 1, ci_lty = 1,
                        xline = NULL, yline = NULL, ref_color = "gray40", ref_lty = 2, ref_lwd = 1,
                        xlab = NULL, ylab = NULL, xtitle = NULL, ytitle = NULL, title = NULL,
                        subtitle = NULL, note = NULL, color = NULL, palette = "journal", alpha = 1,
                        legend = TRUE, legend_position = "topright", theme = "journal", size = 11,
                        file = NULL, width = 7, height = 5, dpi = 300, show = TRUE, bg = "white",
                        vline = NULL, hline = NULL) {
  call <- match.call(); env <- parent.frame()
  refs <- .r4vn_reference_lines(xline, yline, vline, hline)
  obj <- if (missing(result) || is.null(result)) .r4vn_get_last_margins(TRUE) else result
  if (!inherits(obj, "r4vn_margins") || is.null(obj$raw$margins)) stop("`result` must be returned by margins().", call. = FALSE)
  d <- obj$raw$margins
  candidates <- setdiff(names(d), c("Margin", "SE", "Lower", "Upper", "N"))
  x_expr <- substitute(x); by_expr <- substitute(by)
  xname <- if (.r4vn_expr_is_null(x_expr)) {
    varying <- candidates[vapply(d[candidates], function(z) length(unique(z[!is.na(z)])) > 1L, logical(1))]
    if (length(varying)) varying[1L] else if (length(candidates)) candidates[1L] else NULL
  } else .r4vn_resolve_name_spec(x_expr, d, env, "x", multiple = FALSE)
  if (is.null(xname)) { d$.Scenario <- seq_len(nrow(d)); xname <- ".Scenario" }
  byname <- if (.r4vn_expr_is_null(by_expr)) NULL else .r4vn_resolve_name_spec(by_expr, d, env, "by", multiple = FALSE)
  xv <- d[[xname]]; numeric_x <- is.numeric(xv)
  xp <- if (numeric_x) xv else as.numeric(factor(xv, levels = unique(xv)))
  groups <- if (is.null(byname)) factor(rep("Margin", nrow(d))) else factor(d[[byname]], levels = unique(d[[byname]]))
  cols <- .r4vn_colors(nlevels(groups), color, palette, alpha)
  ci_cols <- if (is.null(ci_color)) vapply(cols, grDevices::adjustcolor, character(1), alpha.f = ci_alpha) else .r4vn_colors(nlevels(groups), ci_color, palette, ci_alpha)
  line_type <- rep(line_type, length.out = nlevels(groups))
  point_pch <- rep(point_pch, length.out = nlevels(groups))
  if (!is.numeric(line_width) || length(line_width) != 1L || is.na(line_width) || line_width <= 0) stop("`line_width` must be positive.", call. = FALSE)
  if (!is.numeric(point_size) || length(point_size) != 1L || is.na(point_size) || point_size <= 0) stop("`point_size` must be positive.", call. = FALSE)
  if (!is.numeric(ci_alpha) || length(ci_alpha) != 1L || is.na(ci_alpha) || ci_alpha < 0 || ci_alpha > 1) stop("`ci_alpha` must be between 0 and 1.", call. = FALSE)
  if (!is.numeric(ci_lwd) || length(ci_lwd) != 1L || is.na(ci_lwd) || ci_lwd <= 0) stop("`ci_lwd` must be positive.", call. = FALSE)
  xt <- xtitle %||% xname; yt <- ytitle %||% paste0("Margin (", obj$raw$type, ")")
  draw <- function() {
    old <- .r4vn_theme(theme, size, note); on.exit(graphics::par(old), add = TRUE)
    yr <- range(c(d$Lower, d$Upper, d$Margin), finite = TRUE)
    graphics::plot(range(xp, finite = TRUE), yr, type = "n", axes = FALSE, xlab = xt, ylab = yt, main = title)
    if (numeric_x) .r4vn_axis(1, xp, xlab, NULL) else graphics::axis(1, at = seq_along(unique(xv)), labels = if (is.null(xlab)) unique(xv) else xlab)
    .r4vn_axis(2, d$Margin, ylab, NULL); graphics::box(bty = graphics::par("bty"))
    for (i in seq_len(nlevels(groups))) {
      ii <- which(groups == levels(groups)[i]); oo <- order(xp[ii]); jj <- ii[oo]
      if (isTRUE(ci)) graphics::segments(xp[jj], d$Lower[jj], xp[jj], d$Upper[jj], col = ci_cols[i], lwd = ci_lwd, lty = ci_lty)
      if (isTRUE(line) && length(jj) > 1L) graphics::lines(xp[jj], d$Margin[jj], col = cols[i], lwd = line_width, lty = line_type[i])
      if (isTRUE(points)) graphics::points(xp[jj], d$Margin[jj], col = cols[i], pch = point_pch[i], cex = point_size)
    }
    .r4vn_add_reference_lines(refs$xline, refs$yline, ref_color, ref_lty, ref_lwd)
    if (!is.null(byname) && isTRUE(legend)) graphics::legend(legend_position, legend = levels(groups), col = cols, lwd = line_width, lty = line_type, pch = if (isTRUE(points)) point_pch else NA, pt.cex = point_size, bty = "n")
    .r4vn_add_titles(subtitle, note, size)
  }
  .r4vn_render(draw, file, show, width, height, dpi, bg)
  .r4vn_graph_result("margins", d, call, file, draw = draw)
}

.r4vn_predict_type <- function(type, model) {
  key <- tolower(trimws(as.character(type)[1L]))
  common <- c(
    fitted = "response", predicted = "response", yhat = "response", mu = "response",
    xb = "link", linear = "link", linear.predictor = "link",
    resid = "residual", rawres = "residual",
    standardized = "standardized", standardised = "standardized", stdres = "standardized", rstandard = "standardized",
    studentized = "studentized", studentised = "studentized", studres = "studentized", rstudent = "studentized",
    hat = "leverage", hatvalue = "leverage", hatvalues = "leverage",
    cook = "cooksd", cooks = "cooksd", cooks.d = "cooksd", cooks.distance = "cooksd",
    se = "se.fit", sefit = "se.fit",
    lcl = "lower", ucl = "upper"
  )
  if (key %in% names(common)) key <- common[[key]]
  if (inherits(model, "coxph")) {
    aliases <- c(pr = "risk", probability = "risk", response = "risk", risk = "risk",
                 fitted = "risk", link = "lp", xb = "lp", linear = "lp",
                 linear.predictor = "lp", resid = "martingale")
    if (key %in% names(aliases)) key <- aliases[[key]]
    if (key == "auto") key <- "risk"
    return(key)
  }
  if (key %in% c("pr", "probability")) key <- "response"
  if (inherits(model, "lm") && !inherits(model, "glm") && key == "link") key <- "response"
  if (key == "auto") key <- "response"
  key
}

.r4vn_predict_matrix_column <- function(z, term = NULL, type = "diagnostic") {
  if (!is.matrix(z) && !is.data.frame(z)) return(as.numeric(z))
  z <- as.matrix(z)
  if (ncol(z) == 1L) return(as.numeric(z[, 1L]))
  if (is.null(term)) {
    stop("Prediction type `", type, "` returns multiple columns. Supply `term` to select one. Available terms: ",
         paste(colnames(z) %||% seq_len(ncol(z)), collapse = ", "), ".", call. = FALSE)
  }
  if (is.numeric(term) && length(term) == 1L && is.finite(term)) {
    j <- as.integer(term)
    if (j < 1L || j > ncol(z)) stop("`term` column index is out of range.", call. = FALSE)
    return(as.numeric(z[, j]))
  }
  nm <- as.character(term)[1L]
  if (is.null(colnames(z)) || !nm %in% colnames(z)) {
    stop("`term = ", shQuote(nm), "` was not found. Available terms: ", paste(colnames(z) %||% seq_len(ncol(z)), collapse = ", "), ".", call. = FALSE)
  }
  as.numeric(z[, nm])
}

.r4vn_postestimate_value <- function(fit, ptype, newdata = NULL, term = NULL, dots = list()) {
  # Residual and influence statistics are defined on the estimation sample.
  if (ptype == "standardized") return(.r4vn_predict_matrix_column(tryCatch(stats::rstandard(fit), error = function(e) stop("Standardized residuals are not available for this model: ", conditionMessage(e), call. = FALSE)), term, ptype))
  if (ptype == "studentized") return(.r4vn_predict_matrix_column(tryCatch(stats::rstudent(fit), error = function(e) stop("Studentized residuals are not available for this model: ", conditionMessage(e), call. = FALSE)), term, ptype))
  if (ptype == "leverage") return(.r4vn_predict_matrix_column(tryCatch(stats::hatvalues(fit), error = function(e) stop("Leverage values are not available for this model: ", conditionMessage(e), call. = FALSE)), term, ptype))
  if (ptype == "cooksd") return(.r4vn_predict_matrix_column(tryCatch(stats::cooks.distance(fit), error = function(e) stop("Cook's distance is not available for this model: ", conditionMessage(e), call. = FALSE)), term, ptype))
  if (ptype == "dffits") return(.r4vn_predict_matrix_column(tryCatch(stats::dffits(fit), error = function(e) stop("DFFITS is not available for this model: ", conditionMessage(e), call. = FALSE)), term, ptype))
  if (ptype == "covratio") return(.r4vn_predict_matrix_column(tryCatch(stats::covratio(fit), error = function(e) stop("COVRATIO is not available for this model: ", conditionMessage(e), call. = FALSE)), term, ptype))
  if (ptype == "dfbeta") {
    z <- if (inherits(fit, "coxph")) stats::residuals(fit, type = "dfbeta") else stats::dfbeta(fit)
    return(.r4vn_predict_matrix_column(z, term, ptype))
  }
  if (ptype == "dfbetas") {
    z <- if (inherits(fit, "coxph")) stats::residuals(fit, type = "dfbetas") else stats::dfbetas(fit)
    return(.r4vn_predict_matrix_column(z, term, ptype))
  }

  residual_types <- c("residual", "pearson", "deviance", "working", "martingale", "score", "schoenfeld", "scaledsch", "partial")
  if (ptype %in% residual_types) {
    rt <- if (ptype == "residual") "response" else ptype
    z <- tryCatch(stats::residuals(fit, type = rt), error = function(e) {
      if (ptype == "residual") stats::residuals(fit) else stop("Residual type `", ptype, "` is not available for this model: ", conditionMessage(e), call. = FALSE)
    })
    return(.r4vn_predict_matrix_column(z, term, ptype))
  }

  prediction_data <- if (is.null(newdata)) NULL else .r4vn_prepare_prediction_data(fit, newdata)
  if (ptype == "se.fit") {
    args <- c(list(object = fit), if (is.null(prediction_data)) list() else list(newdata = prediction_data), dots, list(se.fit = TRUE))
    if (inherits(fit, "glm") && is.null(args$type)) args$type <- "response"
    z <- tryCatch(do.call(stats::predict, args), error = function(e) stop("Prediction standard errors are not available: ", conditionMessage(e), call. = FALSE))
    if (!is.list(z) || is.null(z$se.fit)) stop("This model's predict method does not return `se.fit`.", call. = FALSE)
    return(as.numeric(z$se.fit))
  }
  if (ptype %in% c("lower", "upper")) {
    if (!inherits(fit, "lm") || inherits(fit, "glm")) stop("`lower` and `upper` are currently available for linear `lm` models. For GLMs use `se.fit` and model-appropriate transformation.", call. = FALSE)
    if (is.null(dots$interval)) dots$interval <- "confidence"
    args <- c(list(object = fit), if (is.null(prediction_data)) list() else list(newdata = prediction_data), dots)
    z <- do.call(stats::predict, args)
    if (!is.matrix(z)) stop("This prediction did not return interval columns.", call. = FALSE)
    col <- if (ptype == "lower") "lwr" else "upr"
    if (!col %in% colnames(z)) stop("The requested interval column was not returned.", call. = FALSE)
    return(as.numeric(z[, col]))
  }

  args <- c(list(object = fit), if (is.null(prediction_data)) list() else list(newdata = prediction_data), dots)
  if (inherits(fit, "coxph")) {
    args$type <- ptype
    if (ptype %in% c("risk", "lp") && is.null(args$reference)) args$reference <- "zero"
  } else if (inherits(fit, "glm")) {
    args$type <- ptype
  } else if (!inherits(fit, "rq") && ptype != "response") {
    args$type <- ptype
  }
  z <- tryCatch(do.call(stats::predict, args), error = function(e) stop("Prediction failed: ", conditionMessage(e), call. = FALSE))
  .r4vn_predict_matrix_column(z, term, ptype)
}

.r4vn_align_estimation_vector <- function(model, z, n, data) {
  z <- as.numeric(z)
  if (length(z) == n) return(z)
  out <- rep(NA_real_, n)
  rn_model <- tryCatch(rownames(stats::model.frame(model)), error = function(e) NULL)
  rn_data <- rownames(data)
  if (!is.null(rn_model) && !is.null(rn_data) && length(rn_model) == length(z)) {
    hit <- match(rn_model, rn_data); good <- !is.na(hit); out[hit[good]] <- z[good]; return(out)
  }
  used <- tryCatch(.r4vn_used_rows(model, n), error = function(e) NULL)
  if (!is.null(used) && length(used) == length(z)) out[used] <- z
  out
}

.r4vn_write_prediction <- function(d, name, values, data_expr, env, active_mode) {
  d[[name]] <- values
  if (isTRUE(active_mode)) {
    link <- .r4vn_active_link()
    .r4vn_set_active(d, name = .r4vn_active_name(), source = .r4vn_active_source(), quiet = TRUE,
                     object_name = if (is.null(link)) NULL else link$name,
                     object_env = if (is.null(link)) NULL else link$env)
  } else if (is.symbol(data_expr)) {
    nm <- as.character(data_expr)
    target <- .r4vn_find_binding_env(nm, env)
    if (!is.null(target) && !bindingIsLocked(nm, target)) assign(nm, d, envir = target)
  }
  invisible(d)
}

#' Prediction and postestimation diagnostics
#' @usage
#' predict(
#'   object = NULL, ..., newvar = NULL, type = "auto", term = NULL,
#'   data = NULL, replace = FALSE, show = TRUE
#' )
#'
#' @description
#' `predict()` provides one consistent R4VN postestimation interface. With an
#' explicit fitted model and no `newvar`, ordinary prediction types continue to
#' delegate to the model's `stats::predict()` method. R4VN also adds common
#' residual and influence statistics. In variable-generation mode, provide
#' `newvar` and R4VN writes the selected statistic back to `data` or the active
#' data frame while preserving omitted estimation rows as `NA`.
#'
#' @param object Optional fitted model or R4VN model result. If omitted when
#'   generating a variable, the most recent active model is used.
#' @param ... Additional model-specific arguments. `newdata` may be supplied
#'   here for ordinary explicit-model prediction and is normalized to the
#'   predictor types retained by the fitted model. Other examples include
#'   `level` and `interval` for linear-model prediction or arguments accepted
#'   by the underlying model's prediction method.
#' @param newvar Name of a variable to create. It may be unquoted, for example
#'   `newvar = stdres`, or supplied as one character string. If omitted, the
#'   requested statistic is returned instead of being written to data.
#' @param type Statistic to obtain. Common prediction aliases are
#'   `"auto"`, `"response"`, `"fitted"`, `"predicted"`, `"probability"`/`"pr"`,
#'   and `"link"`/`"xb"`. Residual types include `"residual"`, `"pearson"`,
#'   `"deviance"`, `"working"`, `"standardized"`/`"stdres"`/`"rstandard"`, and
#'   `"studentized"`/`"studres"`/`"rstudent"`. Influence statistics include
#'   `"leverage"`/`"hat"`, `"cooksd"`, `"dffits"`, `"covratio"`, `"dfbeta"`,
#'   and `"dfbetas"`. `"se.fit"` returns prediction standard errors. For linear
#'   models, `"lower"` and `"upper"` return confidence-limit columns. Cox models
#'   additionally support `"risk"`, `"lp"`, `"expected"`, `"terms"`,
#'   `"martingale"`, `"deviance"`, `"score"`, `"schoenfeld"`, `"scaledsch"`,
#'   and `"partial"` when supported by `survival`.
#' @param term Optional coefficient/term name or column number when a statistic
#'   naturally returns several columns, notably `dfbeta`, `dfbetas`, Cox score,
#'   Schoenfeld, scaled Schoenfeld, partial residuals, or term predictions.
#'   If omitted for a multi-column result, R4VN reports the available terms.
#' @param data Data frame used for prediction and/or receiving `newvar`. When
#'   omitted in generation mode, the active data frame is used.
#' @param replace Logical; allow an existing `newvar` to be overwritten.
#' @param show Logical; display a short generation message. Default `TRUE`.
#'
#' @details
#' Standardized residuals are computed with `stats::rstandard()` and
#' studentized residuals with `stats::rstudent()` when those methods are
#' available. These are different from raw residuals. For linear regression,
#' leverage is obtained with `hatvalues()`, Cook's distance with
#' `cooks.distance()`, DFFITS with `dffits()`, and COVRATIO with `covratio()`.
#'
#' Influence statistics and residuals are defined for the estimation sample.
#' When they are written to the original active data, observations omitted from
#' model fitting because of missing values are filled with `NA`.
#'
#' When R4VN compact syntax declared a predictor categorical (for example
#' `i.htn`) but the original data store it as numeric 0/1, prediction data are
#' automatically reconstructed with the factor levels retained by the fitted
#' model. Unknown new levels remain an error rather than being silently recoded.
#'
#' @return Without `newvar`, returns the requested prediction, residual, or
#'   diagnostic statistic. With `newvar`, invisibly returns the updated data
#'   frame after writing the generated variable.
#'
#' @examples
#' # Linear regression: fitted values and regression diagnostics
#' d <- data.frame(
#'   y = c(12, 15, 17, 20, 21, 25, 28, 31, 35, 38),
#'   age = seq(20, 65, by = 5),
#'   bmi = c(19, 21, 20, 23, 25, 24, 27, 28, 30, 29)
#' )
#' usedf(d)
#' m1 <- regress(y, c.age, c.bmi, show = FALSE)
#' predict(m1, type = "response")
#' predict(m1, type = "standardized")
#' predict(m1, type = "studentized")
#' predict(m1, type = "leverage")
#' predict(m1, type = "cooksd")
#' predict(m1, type = "dffits")
#' predict(m1, type = "covratio")
#'
#' # Store diagnostics in the active data frame
#' predict(newvar = fitted_y, type = "fitted", show = FALSE)
#' predict(newvar = residual_y, type = "residual", show = FALSE)
#' predict(newvar = stdres, type = "standardized", show = FALSE)
#' predict(newvar = studres, type = "studentized", show = FALSE)
#' predict(newvar = leverage, type = "leverage", show = FALSE)
#' predict(newvar = cooksd, type = "cooksd", show = FALSE)
#'
#' # DFBETA/DFBETAS are coefficient-specific; select a term when generating
#' names(stats::coef(m1$raw$model))
#' predict(m1, type = "dfbetas", term = "age")
#' predict(m1, newvar = dfb_age, type = "dfbetas", term = "age", show = FALSE)
#'
#' # Linear-model prediction standard error and confidence limits
#' predict(m1, type = "se.fit")
#' predict(m1, type = "lower", level = 0.95)
#' predict(m1, type = "upper", level = 0.95)
#'
#' # Logistic regression: probability and residual diagnostics
#' g <- data.frame(
#'   outcome = factor(c(0,0,0,0,1,0,1,1,1,1,1,1), levels = 0:1,
#'                    labels = c("No", "Yes")),
#'   age = seq(25, 80, by = 5),
#'   bmi = c(20,21,22,24,23,26,25,28,29,30,31,33)
#' )
#' usedf(g)
#' m2 <- logistic(outcome, c.age, c.bmi, event = "Yes", show = FALSE)
#' predict(m2, type = "probability")
#' predict(m2, type = "pearson")
#' predict(m2, type = "deviance")
#' predict(m2, type = "standardized")
#' predict(m2, type = "leverage")
#' predict(m2, type = "cooksd")
#' @export
predict <- function(object = NULL, ..., newvar = NULL, type = "auto", term = NULL,
                    data = NULL, replace = FALSE, show = TRUE) {
  env <- parent.frame()
  oexpr <- substitute(object)
  nvexpr <- substitute(newvar)
  dexpr <- substitute(data)
  object_missing <- missing(object) || .r4vn_expr_is_null(oexpr)
  newvar_missing <- missing(newvar) || .r4vn_expr_is_null(nvexpr)
  type_missing <- missing(type) || identical(type, "auto")
  dots <- list(...)

  obj <- NULL
  if (!object_missing) obj <- tryCatch(eval(oexpr, envir = env), error = function(e) NULL)

  # Explicit model and no generated variable: preserve base prediction behavior
  # for the ordinary default, but support all R4VN diagnostic aliases when a
  # type is explicitly requested.
  if (!object_missing && !is.null(obj) && newvar_missing) {
    fit <- .r4vn_unwrap_model(obj)

    # `newdata` is conventionally supplied through ... because the R4VN
    # generation interface reserves `data`. Normalize either spelling here so
    # compact factor declarations such as i.htn are reconstructed before the
    # underlying stats::predict() method sees the data.
    if (!missing(data) && !is.null(data) && !is.null(dots$newdata)) {
      stop("Supply only one of `data` or `newdata` for prediction, not both.", call. = FALSE)
    }
    if (!missing(data) && !is.null(data)) {
      dots$newdata <- .r4vn_prepare_prediction_data(fit, data)
    } else if (!is.null(dots$newdata)) {
      dots$newdata <- .r4vn_prepare_prediction_data(fit, dots$newdata)
    }

    if (type_missing) {
      args <- c(list(object = fit), dots)
      return(do.call(stats::predict, args))
    }
    ptype <- .r4vn_predict_type(type, fit)
    special <- c("standardized", "studentized", "leverage", "cooksd", "dffits", "covratio", "dfbeta", "dfbetas",
                 "residual", "pearson", "deviance", "working", "martingale", "score", "schoenfeld", "scaledsch", "partial",
                 "se.fit", "lower", "upper")
    if (ptype %in% special) {
      nd <- if (!is.null(dots$newdata)) dots$newdata else NULL
      dots$newdata <- NULL
      return(.r4vn_postestimate_value(fit, ptype, newdata = nd, term = term, dots = dots))
    }
    args <- c(list(object = fit), dots)
    if (inherits(fit, "coxph")) {
      args$type <- ptype
      if (ptype %in% c("risk", "lp") && is.null(args$reference)) args$reference <- "zero"
    } else if (inherits(fit, "glm")) args$type <- ptype
    else if (!inherits(fit, "rq") && ptype != "response") args$type <- ptype
    return(do.call(stats::predict, args))
  }

  # Stata-like shorthand: predict(phat, type = "pr") when phat is not an
  # existing object means generate phat from the active model.
  if (!object_missing && is.null(obj) && newvar_missing && is.symbol(oexpr)) {
    varname <- as.character(oexpr)
    fit <- .r4vn_get_active_model(TRUE)
  } else {
    fit <- if (!object_missing && !is.null(obj)) .r4vn_unwrap_model(obj) else .r4vn_get_active_model(TRUE)
    if (newvar_missing) stop("Supply `newvar`, for example `predict(newvar = phat, type = 'probability')`.", call. = FALSE)
    varname <- if (is.symbol(nvexpr)) as.character(nvexpr) else {
      vv <- eval(nvexpr, envir = env)
      if (!is.character(vv) || length(vv) != 1L || !nzchar(vv)) stop("`newvar` must be a name or one non-empty character string.", call. = FALSE)
      vv
    }
  }

  active_mode <- missing(data) || is.null(data)
  d <- if (active_mode) .r4vn_get_active(TRUE) else data
  if (!is.data.frame(d)) stop("`data` must be a data frame.", call. = FALSE)
  if (varname %in% names(d) && !isTRUE(replace)) stop("Variable `", varname, "` already exists; use `replace = TRUE` to overwrite it.", call. = FALSE)

  ptype <- .r4vn_predict_type(type, fit)
  estimation_types <- c("standardized", "studentized", "leverage", "cooksd", "dffits", "covratio", "dfbeta", "dfbetas",
                        "residual", "pearson", "deviance", "working", "martingale", "score", "schoenfeld", "scaledsch", "partial")
  if (ptype %in% estimation_types) {
    z <- .r4vn_postestimate_value(fit, ptype, newdata = NULL, term = term, dots = dots)
    z <- .r4vn_align_estimation_vector(fit, z, nrow(d), d)
  } else {
    z <- .r4vn_postestimate_value(fit, ptype, newdata = d, term = term, dots = dots)
    if (length(z) != nrow(d)) z <- .r4vn_align_estimation_vector(fit, z, nrow(d), d)
  }

  out <- .r4vn_write_prediction(d, varname, z, dexpr, env, active_mode)
  if (isTRUE(show)) message("Generated `", varname, "` from the active/fitted model (type = ", ptype, ").")
  out
}

.r4vn_lincom_vector <- function(expr, coef_names, env) {
  if (is.character(expr) && length(expr) == 1L) expr <- parse(text = expr, keep.source = FALSE)[[1L]]
  p <- length(coef_names)
  evaluate_at <- function(values) {
    ee <- new.env(parent = env)
    for (i in seq_along(coef_names)) assign(coef_names[i], values[i], envir = ee)
    z <- eval(expr, envir = ee)
    if (!is.numeric(z) || length(z) != 1L || !is.finite(z)) {
      stop("`lincom` expression must evaluate to one finite numeric linear combination of coefficient names.", call. = FALSE)
    }
    as.numeric(z)
  }
  zero <- numeric(p)
  constant <- evaluate_at(zero)
  L <- numeric(p)
  for (i in seq_len(p)) {
    one <- zero; one[i] <- 1
    L[i] <- evaluate_at(one) - constant
  }
  # Guard against nonlinear expressions (e.g. b1*b2 or b1^2). A deterministic
  # nontrivial point is enough to verify that the affine reconstruction holds.
  probe <- seq_len(p) / (p + 1)
  actual <- evaluate_at(probe)
  reconstructed <- constant + sum(L * probe)
  tol <- 1e-8 * max(1, abs(actual), abs(reconstructed))
  if (!isTRUE(all.equal(actual, reconstructed, tolerance = tol, check.attributes = FALSE))) {
    stop("`lincom` supports linear/affine combinations only; nonlinear coefficient expressions are not allowed.", call. = FALSE)
  }
  list(L = L, constant = constant, expression = expr)
}

#' Linear combinations of fitted-model coefficients
#' @usage
#' lincom(
#'   ..., model = NULL, rhs = 0, exp = FALSE, level = 0.95, digits = 3,
#'   p_digits = 3, show = TRUE, console = FALSE
#' )
#'
#' @param ... One or more unquoted expressions or character expressions using
#'   coefficient names, e.g. `lincom(age + 2*bmi)`.
#' @param model Optional model; the active model is used by default.
#' @param rhs Null value for the linear combination.
#' @param exp Exponentiate the estimate and confidence interval.
#' @param level Confidence level.
#' @param digits,p_digits Formatting controls.
#' @param show,console Display controls.
#' @return An `r4vn_stat` object.
#' @examples
#' d <- data.frame(
#'   y = c(50, 54, 57, 61, 65, 68, 72, 76),
#'   age = c(20, 25, 30, 35, 40, 45, 50, 55),
#'   bmi = c(20, 22, 21, 24, 25, 27, 26, 29)
#' )
#' regress(y, c.age, c.bmi, data = d, show = FALSE)
#' lincom(age + 2 * bmi, show = FALSE)
#' lincom("age - bmi", show = FALSE)
#' @export
lincom <- function(..., model = NULL, rhs = 0, exp = FALSE, level = 0.95,
                   digits = 3, p_digits = 3, show = TRUE, console = FALSE) {
  call <- match.call(); env <- parent.frame(); exprs <- as.list(substitute(list(...)))[-1L]
  if (!length(exprs)) stop("Supply at least one linear-combination expression.", call. = FALSE)
  fit <- if (missing(model) || is.null(model)) .r4vn_get_active_model(TRUE) else .r4vn_unwrap_model(model)
  b <- stats::coef(fit); V <- .r4vn_post_vcov(fit)
  if (is.null(V)) stop("A covariance matrix is not available for this model.", call. = FALSE)
  good <- is.finite(b) & names(b) %in% rownames(V) & names(b) %in% colnames(V)
  b <- b[good]; V <- V[names(b), names(b), drop = FALSE]
  use_t <- inherits(fit, "lm") && !inherits(fit, "glm") && is.finite(stats::df.residual(fit)) && stats::df.residual(fit) > 0
  crit <- if (use_t) stats::qt(1 - (1 - level) / 2, df = stats::df.residual(fit)) else stats::qnorm(1 - (1 - level) / 2)
  rows <- vector("list", length(exprs))
  for (i in seq_along(exprs)) {
    ex <- exprs[[i]]
    # Character linear combinations may be supplied directly, through a
    # character object, or through an expression that constructs one (for
    # example paste0("`", coef_name, "` + 1")).  Force the caller-side
    # expression once and only replace `ex` when it yields a scalar string.
    # Numeric results are deliberately ignored here so ordinary unquoted
    # coefficient expressions such as age + 2*bmi are still evaluated in the
    # coefficient environment below rather than in the caller's data scope.
    val <- tryCatch(eval(ex, env), error = function(e) NULL)
    if (is.character(val) && length(val) == 1L && !is.na(val)) ex <- val
    lc <- .r4vn_lincom_vector(ex, names(b), env)
    L <- lc$L
    est <- lc$constant + sum(L * b)
    vv <- drop(t(L) %*% V %*% L)
    se <- if (is.finite(vv) && vv >= 0) sqrt(vv) else NA_real_
    stat <- if (is.finite(se) && se > 0) (est - rhs) / se else if (isTRUE(all.equal(est, rhs))) 0 else sign(est - rhs) * Inf
    p <- if (is.finite(stat)) {
      if (use_t) 2 * stats::pt(abs(stat), df = stats::df.residual(fit), lower.tail = FALSE) else 2 * stats::pnorm(abs(stat), lower.tail = FALSE)
    } else if (is.infinite(stat)) 0 else NA_real_
    ci <- if (is.finite(se)) est + c(-1, 1) * crit * se else c(NA_real_, NA_real_)
    disp <- if (isTRUE(exp)) c(exp(est), exp(ci)) else c(est, ci)
    se_disp <- if (isTRUE(exp) && is.finite(se)) exp(est) * se else se
    rows[[i]] <- data.frame(Expression = paste(deparse(ex, width.cutoff = 500L), collapse = ""),
      Estimate = .r4vn_num(disp[1L], digits), SE = .r4vn_num(se_disp, digits),
      Lower = .r4vn_num(disp[2L], digits), Upper = .r4vn_num(disp[3L], digits),
      Statistic = .r4vn_num(stat, digits), p = .r4vn_p(p, p_digits), stringsAsFactors = FALSE, check.names = FALSE)
  }
  tab <- do.call(rbind, rows); rownames(tab) <- NULL
  .r4vn_show(.r4vn_result("Linear combination of model coefficients", list("Linear combinations" = tab),
    raw = list(model = fit, table = tab, vcov = V), call = call), show = show, console = console)
}