R/advanced-models.R

In R4VN: Health Data Analysis and Publication-Ready Reporting

# ==========================================================================
# Advanced regression and nonparametric trend commands
# ==========================================================================

#' Quantile regression
#' @usage qregress(y, ..., vars = NULL, data = NULL, tau = 0.5, method = "br", se = "nid", weights = NULL, subset = NULL, ref = NULL, noconstant = FALSE, diagnosis = FALSE, level = 0.95, digits = 3, p_digits = 3, show = TRUE, console = FALSE)
#'
#' @description
#' Fits one or more conditional quantiles using the `quantreg` package while
#' retaining the R4VN model syntax (`vars()`, factor directives, active data,
#' reference levels, weights, and publication-ready tables).
#'
#' @param y Outcome variable or regression formula.
#' @param ... Predictor terms.
#' @param vars Optional R4VN `vars(...)` predictor specification.
#' @param data Data frame; active R4VN data is used when omitted.
#' @param tau Quantile(s) between 0 and 1, e.g. `.5` or `c(.25,.5,.75)`.
#' @param method Algorithm passed to `quantreg::rq()`.
#' @param se Standard-error method passed to `summary.rq()`; common choices are
#'   `"nid"`, `"iid"`, `"ker"`, `"boot"`, and `"rank"`.
#' @param weights,subset,ref Model controls consistent with other R4VN models.
#' @param noconstant Remove the intercept.
#' @param diagnosis Logical; if `TRUE`, append quantile-regression diagnostics for the primary fitted quantile, including residual median/MAD/IQR, residual balance around zero, and quantile check loss. Default `FALSE`.
#' @param level Confidence level.
#' @param digits,p_digits Formatting controls.
#' @param show,console Display controls.
#' @return An `r4vn_stat`; `raw$model` is the primary fitted `rq` model and
#'   `raw$models` contains all requested quantiles.
#' @export
#' @examples
#' if (requireNamespace("quantreg", quietly = TRUE)) {
#'   # Use a reasonably sized, full-rank data set so the example is stable
#'   # across quantreg and R versions.
#'   d <- datasets::mtcars
#'   d$am <- factor(d$am, levels = c(0, 1),
#'                  labels = c("Automatic", "Manual"))
#'
#'   # Median regression with one continuous and one categorical predictor.
#'   qregress(mpg, c.wt, i.am, data = d, tau = .5,
#'            se = "iid", show = FALSE)
#'
#'   # Fit several conditional quantiles in one call.
#'   qregress(mpg, c.wt, i.am, data = d,
#'            tau = c(.25, .5, .75), se = "iid", show = FALSE)
#'
#'   # Request the R4VN quantile-regression diagnostic section.
#'   qregress(mpg, c.wt, i.am, data = d, tau = .5,
#'            se = "iid", diagnosis = TRUE, show = FALSE)
#' }
qregress <- function(y, ..., vars = NULL, data = NULL, tau = 0.5,
                     method = "br", se = "nid", weights = NULL, subset = NULL,
                     ref = NULL, noconstant = FALSE, diagnosis = FALSE, level = 0.95,
                     digits = 3, p_digits = 3, show = TRUE, console = FALSE) {
  if (!requireNamespace("quantreg", quietly = TRUE)) {
    stop("Package `quantreg` is required for qregress(). Install it with install.packages('quantreg').", call. = FALSE)
  }
  if (!is.numeric(tau) || !length(tau) || any(!is.finite(tau) | tau <= 0 | tau >= 1)) stop("`tau` must contain values strictly between 0 and 1.", call. = FALSE)
  if (!is.numeric(level) || length(level) != 1L || level <= 0 || level >= 1) stop("`level` must be between 0 and 1.", call. = FALSE)
  env <- parent.frame(); 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), 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)

  fits <- vector("list", length(tau)); names(fits) <- format(tau, trim = TRUE)
  sections <- list(); raw_coef <- list(); covs <- list(); zcrit <- stats::qnorm(1 - (1 - level) / 2)
  for (i in seq_along(tau)) {
    args <- list(formula = fit_formula, tau = tau[i], data = prep$data, method = method, na.action = stats::na.omit)
    if (!is.null(prep$weights)) args$weights <- prep$weights
    fit <- do.call(quantreg::rq, args); fits[[i]] <- fit
    sm <- tryCatch(summary(fit, se = se, covariance = TRUE), error = function(e) summary(fit, se = se))
    cm <- as.matrix(sm$coefficients)
    est <- cm[, 1L]
    cnm <- tolower(colnames(cm) %||% rep("", ncol(cm)))
    # Most summary.rq() SE methods return Value/Std. Error/t value/Pr(>|t|).
    # Rank inversion instead returns Value/Lower Bd/Upper Bd; preserve those
    # confidence limits rather than misreading them as an SE and test statistic.
    lower_col <- grep("lower", cnm)[1L]
    upper_col <- grep("upper", cnm)[1L]
    rank_bounds <- length(lower_col) && length(upper_col) && is.finite(lower_col) && is.finite(upper_col)
    if (isTRUE(rank_bounds)) {
      sev <- rep(NA_real_, length(est))
      stat <- rep(NA_real_, length(est))
      pp <- rep(NA_real_, length(est))
      ci <- cbind(cm[, lower_col], cm[, upper_col])
    } else {
      se_col <- grep("std|standard.*error", cnm)[1L]
      sev <- if (length(se_col) && is.finite(se_col)) cm[, se_col] else if (ncol(cm) >= 2L) cm[, 2L] else rep(NA_real_, length(est))
      stat_col <- grep("t value|z value|stat", cnm)[1L]
      stat <- if (length(stat_col) && is.finite(stat_col)) cm[, stat_col] else est / sev
      p_col <- grep("pr\\(|p.value|p-value|p value", cnm)[1L]
      pp <- if (length(p_col) && is.finite(p_col)) cm[, p_col] else 2 * stats::pnorm(abs(stat), lower.tail = FALSE)
      ci <- cbind(est - zcrit * sev, est + zcrit * sev)
    }
    tab <- data.frame(Tau = tau[i], Term = sub("^\\(Intercept\\)$", "_cons", names(est)),
                      Coefficient = .r4vn_num(est, digits), SE = .r4vn_num(sev, digits),
                      Lower = .r4vn_num(ci[, 1L], digits), Upper = .r4vn_num(ci[, 2L], digits),
                      Statistic = .r4vn_num(stat, digits), p = .r4vn_p(pp, p_digits),
                      stringsAsFactors = FALSE, check.names = FALSE)
    sections[[paste0("Coefficients: tau = ", format(tau[i], trim = TRUE))]] <- tab
    raw_coef[[i]] <- list(estimate = est, se = sev, statistic = stat, p.value = pp, conf.int = ci)
    covs[[i]] <- if (!is.null(sm$cov)) sm$cov else tryCatch(stats::vcov(fit), error = function(e) NULL)
  }
  primary_i <- which.min(abs(tau - 0.5)); primary <- fits[[primary_i]]
  # quantreg::rq objects do not consistently provide an nobs() method across
  # quantreg/R versions. quantreg itself uses length(residuals(x)) when
  # reporting the number of observations for an rq fit, so use the same
  # portable definition here.
  n_primary <- length(stats::residuals(primary))
  info <- data.frame(Statistic = c("Dependent variable", "Number of observations", "Quantile(s)", "Method", "SE method"),
                     Value = c(.r4vn_deparse1(f[[2L]]), n_primary, paste(format(tau, trim = TRUE), collapse = ", "), method, se),
                     stringsAsFactors = FALSE)
  sections <- c(list("Model summary" = info), sections)
  diagnostics <- NULL
  if (isTRUE(diagnosis)) {
    diagnostics <- .r4vn_model_diagnosis(primary, kind = "quantile", tau = tau[primary_i], digits = digits, p_digits = p_digits)
    sections <- c(sections, diagnostics)
  }
  .r4vn_show(.r4vn_result("Quantile regression", sections,
    raw = list(model = primary, models = fits, vcov = covs[[primary_i]], coefficients = raw_coef,
               tau = tau, method = method, se = se, model.terms = .r4vn_mx_term_labels(primary), diagnostics = diagnostics),
    call = match.call()), show = show, console = console)
}

.r4vn_quote_name <- function(x) paste0("`", gsub("`", "\\`", x, fixed = TRUE), "`")

#' Flexible nonlinear-shape regression using splines or polynomials
#' @usage nlregress(y, x, covariates = NULL, data = NULL, spline = c("natural", "bspline", "polynomial", "linear"), df = 4, degree = 3, knots = NULL, boundary_knots = NULL, family = c("gaussian", "binomial", "poisson"), event = NULL, robust = FALSE, diagnosis = FALSE, level = 0.95, digits = 3, p_digits = 3, show = TRUE, console = FALSE)
#'
#' @description
#' `nlregress()` provides a simple R4VN interface for nonlinear predictor shapes
#' without requiring users to hand-code spline bases. It supports natural cubic
#' splines, B-splines, raw polynomials, or a linear term and can fit Gaussian,
#' binomial, or Poisson outcomes.
#'
#' @param y Outcome variable.
#' @param x Numeric predictor whose functional form is modeled flexibly.
#' @param covariates Optional additional predictors selected with `vars(...)` or
#'   a character vector.
#' @param data Data frame; active R4VN data is used when omitted.
#' @param spline `"natural"`, `"bspline"`, `"polynomial"`, or `"linear"`.
#' @param df Degrees of freedom for spline bases when `knots` is not supplied.
#' @param degree B-spline or polynomial degree.
#' @param knots Optional internal knots on the x scale.
#' @param boundary_knots Optional two boundary knots.
#' @param family `"gaussian"`, `"binomial"`, or `"poisson"`.
#' @param event Event category for a binary binomial or binary Poisson outcome.
#' @param robust Use a sandwich HC0 covariance estimate when available through
#'   R4VN's internal covariance engine.
#' @param diagnosis Logical; if `TRUE`, append model diagnostics appropriate to the selected family: linear diagnostics for Gaussian models, logistic diagnostics for binomial models, and dispersion/goodness-of-fit diagnostics for Poisson models. Default `FALSE`.
#' @param level,digits,p_digits,show,console Standard R4VN controls.
#' @return An `r4vn_stat` with the fitted model in `raw$model`; it therefore
#'   works immediately with [margins()], [predict()], and [lincom()].
#' @examples
#' d <- data.frame(
#'   age = seq(20, 75, by = 5),
#'   bmi = c(20, 21, 22, 24, 23, 25, 26, 27, 29, 28, 30, 31),
#'   sex = factor(rep(c("Female", "Male"), 6)),
#'   y = c(48, 52, 55, 61, 60, 66, 69, 73, 78, 80, 85, 89),
#'   outcome = c(0, 0, 0, 1, 0, 1, 0, 1, 1, 0, 1, 1)
#' )
#' nlregress(y, age, data = d, spline = "natural", df = 3, show = FALSE)
#' nlregress(y, age, covariates = vars(sex, bmi), data = d,
#'           spline = "natural", knots = c(35, 50), show = FALSE)
#' nlregress(outcome, age, data = d, family = "binomial", event = 1,
#'           show = FALSE)
#' nlregress(y, age, data = d, spline = "bspline", df = 4, diagnosis = TRUE, show = FALSE)
#' nlregress(y, age, data = d, spline = "polynomial", degree = 2, show = FALSE)
#' nlregress(outcome, age, data = d, family = "binomial", event = 1, diagnosis = TRUE, show = FALSE)
#' @importFrom splines ns bs
#' @export
nlregress <- function(y, x, covariates = NULL, data = NULL,
                      spline = c("natural", "bspline", "polynomial", "linear"),
                      df = 4, degree = 3, knots = NULL, boundary_knots = NULL,
                      family = c("gaussian", "binomial", "poisson"), event = NULL,
                      robust = FALSE, diagnosis = FALSE, level = 0.95, digits = 3, p_digits = 3,
                      show = TRUE, console = FALSE) {
  call <- match.call(); env <- parent.frame(); d <- .r4vn_stat_data(data)
  spline <- match.arg(spline); family <- match.arg(family)
  if (!is.numeric(df) || length(df) != 1L || !is.finite(df) || df < 1) stop("`df` must be one positive number.", call. = FALSE)
  if (!is.numeric(degree) || length(degree) != 1L || !is.finite(degree) || degree < 1) stop("`degree` must be one positive integer.", call. = FALSE)
  if (!is.null(knots) && (!is.numeric(knots) || any(!is.finite(knots)))) stop("`knots` must contain finite numeric values.", call. = FALSE)
  if (!is.null(boundary_knots) && (!is.numeric(boundary_knots) || length(boundary_knots) != 2L || any(!is.finite(boundary_knots)) || boundary_knots[1L] >= boundary_knots[2L])) stop("`boundary_knots` must contain two increasing finite numeric values.", call. = FALSE)
  yn <- .r4vn_resolve_name_spec(substitute(y), d, env, "y", multiple = FALSE)
  xn <- .r4vn_resolve_name_spec(substitute(x), d, env, "x", multiple = FALSE)
  if (!is.numeric(d[[xn]])) stop("`x` must be numeric.", call. = FALSE)
  cn <- if (missing(covariates) || .r4vn_expr_is_null(substitute(covariates))) character() else .r4vn_resolve_name_spec(substitute(covariates), d, env, "covariates", allow_null = TRUE, multiple = TRUE)
  cn <- setdiff(cn, c(yn, xn))

  xq <- .r4vn_quote_name(xn); yq <- .r4vn_quote_name(yn)
  term <- switch(spline,
    natural = {
      if (is.null(knots)) sprintf("splines::ns(%s, df = %s%s)", xq, as.integer(df), if (is.null(boundary_knots)) "" else paste0(", Boundary.knots = c(", paste(boundary_knots, collapse = ","), ")"))
      else sprintf("splines::ns(%s, knots = c(%s)%s)", xq, paste(knots, collapse = ","), if (is.null(boundary_knots)) "" else paste0(", Boundary.knots = c(", paste(boundary_knots, collapse = ","), ")"))
    },
    bspline = {
      if (is.null(knots)) sprintf("splines::bs(%s, df = %s, degree = %s%s)", xq, as.integer(df), as.integer(degree), if (is.null(boundary_knots)) "" else paste0(", Boundary.knots = c(", paste(boundary_knots, collapse = ","), ")"))
      else sprintf("splines::bs(%s, knots = c(%s), degree = %s%s)", xq, paste(knots, collapse = ","), as.integer(degree), if (is.null(boundary_knots)) "" else paste0(", Boundary.knots = c(", paste(boundary_knots, collapse = ","), ")"))
    },
    polynomial = sprintf("stats::poly(%s, degree = %s, raw = TRUE)", xq, as.integer(degree)),
    linear = xq
  )
  rhs <- c(term, vapply(cn, .r4vn_quote_name, character(1)))
  f <- stats::as.formula(paste(yq, "~", paste(rhs, collapse = " + ")), env = env)

  event_label <- NULL
  if (family %in% c("binomial", "poisson")) {
    vv <- d[[yn]]; lev <- unique(as.character(vv[!is.na(vv)]))
    if (family == "binomial" || !is.null(event) || (!is.numeric(vv) && length(lev) == 2L)) {
      if (length(lev) != 2L) stop("A binary outcome must have exactly two observed values.", call. = FALSE)
      default_event <- if (is.factor(vv)) {
        levels(droplevels(vv))[2L]
      } else if (is.logical(vv)) {
        "TRUE"
      } else if (is.numeric(vv) && all(lev %in% c("0", "1"))) {
        "1"
      } else {
        tail(lev, 1L)
      }
      event_label <- as.character(event %||% default_event)
      if (!event_label %in% lev) stop("`event` was not found in the outcome.", call. = FALSE)
      d[[yn]] <- as.integer(as.character(vv) == event_label)
    }
  }
  fit <- if (family == "gaussian") stats::lm(f, data = d, na.action = stats::na.omit, model = TRUE, x = TRUE, y = TRUE)
         else stats::glm(f, data = d, family = if (family == "binomial") stats::binomial() else stats::poisson(), na.action = stats::na.omit, model = TRUE, x = TRUE, y = TRUE)
  # model.frame() stores the evaluated spline basis rather than the original x.
  # Retain raw predictor columns from the estimation rows so margins(at=...) can
  # vary the original predictor and rebuild the spline basis correctly.
  used_rows <- rownames(stats::model.frame(fit))
  hit_rows <- match(used_rows, rownames(d))
  hit_rows <- hit_rows[!is.na(hit_rows)]
  fit$.r4vn_prediction_data <- d[hit_rows, unique(c(xn, cn)), drop = FALSE]
  V <- if (isTRUE(robust)) tryCatch(.r4vn_model_vcov(fit, "robust", NULL), error = function(e) stats::vcov(fit)) else stats::vcov(fit)
  dist <- if (inherits(fit, "lm") && !inherits(fit, "glm")) "t" else "z"
  cr <- .r4vn_coef_raw(fit, V, level, dist)
  coefs <- .r4vn_coef_table(cr, digits, p_digits, dist, FALSE, "Coefficient")
  info <- data.frame(Statistic = c("Dependent variable", "Flexible predictor", "Functional form", "Number of obs", "AIC", "Event", "Covariance"),
    Value = c(yn, xn, spline, stats::nobs(fit), .r4vn_num(stats::AIC(fit), digits), event_label %||% "", if (robust) "robust" else "model"), stringsAsFactors = FALSE)
  if (family == "gaussian") {
    sm <- summary(fit)
    info <- rbind(info, data.frame(Statistic = c("R-squared", "Adjusted R-squared"), Value = c(.r4vn_num(sm$r.squared, digits), .r4vn_num(sm$adj.r.squared, digits)), stringsAsFactors = FALSE))
  }
  note <- paste0("Flexible term for `", xn, "`: ", term, ".", if (!is.null(knots)) paste0(" Internal knots: ", paste(knots, collapse = ", "), ".") else "")
  sections <- list("Model summary" = info, "Coefficients" = coefs)
  diagnostics <- NULL
  if (isTRUE(diagnosis)) {
    kind <- if (family == "gaussian") "linear" else if (family == "binomial") "logistic" else "poisson"
    diagnostics <- .r4vn_model_diagnosis(fit, kind = kind, digits = digits, p_digits = p_digits)
    sections <- c(sections, diagnostics)
  }
  .r4vn_show(.r4vn_result("Flexible nonlinear-shape regression", sections,
    notes = note, raw = list(model = fit, vcov = V, coefficients = cr, spline = spline, x = xn, knots = knots, boundary.knots = boundary_knots, event = event_label, diagnostics = diagnostics), call = call),
    show = show, console = console)
}

.r4vn_nptrend_cuzick <- function(x, g, scores = NULL) {
  ok <- is.finite(x) & !is.na(g); x <- x[ok]; g <- droplevels(factor(g[ok], ordered = TRUE))
  k <- nlevels(g); if (k < 2L) stop("At least two ordered groups are required.", call. = FALSE)
  sclev <- if (is.null(scores)) seq_len(k) else as.numeric(scores)
  if (length(sclev) != k || any(!is.finite(sclev))) stop("`scores` must contain one finite score per observed group.", call. = FALSE)
  sc <- sclev[as.integer(g)]; r <- rank(x, ties.method = "average")
  Tc <- sum(sc * r); Er <- mean(r); ET <- sum(sc) * Er
  ss_sc <- sum((sc - mean(sc))^2); ss_r <- sum((r - mean(r))^2)
  varT <- ss_sc * ss_r / (length(r) - 1)
  z <- if (varT > 0) (Tc - ET) / sqrt(varT) else NA_real_
  p <- if (is.finite(z)) 2 * stats::pnorm(abs(z), lower.tail = FALSE) else NA_real_
  list(statistic = z, p.value = p, scores = sclev, n = length(r), direction = sign(z), raw = Tc)
}

.r4vn_nptrend_binary <- function(x, g, event = NULL, scores = NULL) {
  ok <- !is.na(x) & !is.na(g); x <- x[ok]; g <- droplevels(factor(g[ok], ordered = TRUE))
  levx <- unique(as.character(x)); if (length(levx) != 2L) stop("Binary trend analysis requires exactly two outcome values.", call. = FALSE)
  ev <- as.character(event %||% tail(levx, 1L)); if (!ev %in% levx) stop("`event` was not found.", call. = FALSE)
  sc <- if (is.null(scores)) seq_len(nlevels(g)) else as.numeric(scores)
  if (length(sc) != nlevels(g)) stop("`scores` must contain one value per group.", call. = FALSE)
  n <- as.numeric(table(g)); e <- as.numeric(tapply(as.character(x) == ev, g, sum))
  fit <- stats::prop.trend.test(e, n, score = sc)
  props <- e / n
  direction <- sign(stats::cor(sc, props, method = "pearson"))
  z <- direction * sqrt(unname(fit$statistic))
  list(statistic = z, chi.square = unname(fit$statistic), p.value = fit$p.value, scores = sc, event = ev, events = e, totals = n, proportions = props)
}

#' Nonparametric test for trend across ordered groups
#' @usage nptrend(x, by, data = NULL, method = c("auto", "cuzick", "cochran-armitage", "spearman", "linear"), event = NULL, scores = NULL, digits = 3, p_digits = 3, show = TRUE, console = FALSE)
#'
#' @description
#' Provides a Cuzick-style rank trend test for quantitative outcomes and a
#' Cochran-Armitage trend test for binary outcomes. Hierarchical R4VN grouping is
#' supported, e.g. `by = vars(province, sex, dose_group)` analyzes the ordered
#' dose trend within province and sex strata.
#'
#' @param x Outcome variable.
#' @param by Ordered grouping variable; with `vars(...)`, the final variable is
#'   the ordered group and preceding variables are strata.
#' @param data Data frame; active data is used when omitted.
#' @param method `"auto"`, `"cuzick"`, `"cochran-armitage"`, `"spearman"`, or
#'   `"linear"`.
#' @param event Event value for a binary outcome.
#' @param scores Optional numeric scores for ordered group levels.
#' @param digits,p_digits,show,console Standard R4VN controls.
#' @return An `r4vn_stat` object.
#' @export
nptrend <- function(x, by, data = NULL,
                    method = c("auto", "cuzick", "cochran-armitage", "spearman", "linear"),
                    event = NULL, scores = NULL, digits = 3, p_digits = 3,
                    show = TRUE, console = FALSE) {
  call <- match.call(); env <- parent.frame(); d <- .r4vn_stat_data(data); method <- match.arg(method)
  xn <- .r4vn_resolve_name_spec(substitute(x), d, env, "x", multiple = FALSE)
  bs <- .r4vn_by_spec(substitute(by), d, env, allow_null = FALSE)
  ids <- .r4vn_strata_indices(d, bs$strata); if (!length(bs$strata)) ids <- list(Overall = seq_len(nrow(d)))
  rows <- list(); raw <- list()
  for (idx in ids) {
    xx <- d[[xn]][idx]; gg0 <- d[[bs$by]][idx]
    gg <- if (is.factor(gg0)) droplevels(gg0) else factor(gg0, levels = unique(gg0[!is.na(gg0)]), ordered = TRUE)
    obs <- xx[!is.na(xx)]; binary <- length(unique(as.character(obs))) == 2L
    binary_hint <- binary && (!is.numeric(xx) || !is.null(event) || all(unique(as.character(obs)) %in% c("0", "1")))
    use <- if (method == "auto") if (binary_hint) "cochran-armitage" else "cuzick" else method
    if (use == "cochran-armitage") z <- .r4vn_nptrend_binary(xx, gg, event, scores)
    else if (use == "cuzick") {
      if (!is.numeric(xx)) stop("Cuzick trend analysis requires a numeric outcome; use `event`/`method = 'cochran-armitage'` for a binary categorical outcome.", call. = FALSE)
      z <- .r4vn_nptrend_cuzick(xx, gg, scores)
    } else {
      if (!is.numeric(xx)) stop("`method = 'spearman'` and `method = 'linear'` require a numeric outcome.", call. = FALSE)
      ok <- !is.na(xx) & !is.na(gg); sclev <- if (is.null(scores)) seq_len(nlevels(droplevels(factor(gg[ok])))) else scores
      gfac <- droplevels(factor(gg[ok])); sc <- sclev[as.integer(gfac)]
      if (use == "spearman") {
        ct <- suppressWarnings(stats::cor.test(as.numeric(xx[ok]), sc, method = "spearman", exact = FALSE)); z <- list(statistic = unname(ct$estimate), p.value = ct$p.value, scores = sclev)
      } else {
        lmfit <- stats::lm(as.numeric(xx[ok]) ~ sc); sm <- summary(lmfit)$coefficients[2L, ]; z <- list(statistic = unname(sm["t value"]), p.value = unname(sm["Pr(>|t|)"]), scores = sclev, slope = unname(sm["Estimate"]))
      }
    }
    lab <- .r4vn_stratum_label(d, bs$strata, idx)
    rows[[length(rows) + 1L]] <- data.frame(Stratum = lab, Method = use,
      Statistic = .r4vn_num(z$statistic, digits), p = .r4vn_p(z$p.value, p_digits),
      Direction = if (is.finite(z$statistic)) if (z$statistic > 0) "Increasing" else if (z$statistic < 0) "Decreasing" else "No direction" else "",
      stringsAsFactors = FALSE, check.names = FALSE)
    raw[[lab]] <- z
  }
  tab <- do.call(rbind, rows); rownames(tab) <- NULL
  .r4vn_show(.r4vn_result("Nonparametric trend test", list("Trend test" = tab),
    notes = paste0("Ordered grouping variable: ", bs$by, ". Group order follows factor level/observed order unless `scores` is supplied."),
    raw = list(results = raw, by = bs, outcome = xn), call = call), show = show, console = console)
}