R/tabscale.R

In R4VN: Health Data Analysis and Publication-Ready Reporting

# =============================================================================
# R4VN::tabscale()
# Comprehensive scale analysis with minimal external dependencies.
#
# Base/recommended R only:
#   - item descriptive statistics
#   - corrected item-total correlations
#   - Cronbach alpha and standardized alpha
#   - alpha if item deleted
#   - split-half/Spearman-Brown, Guttman lambda 6, omega total
#   - KMO, Bartlett test, parallel analysis
#   - principal-axis or maximum-likelihood EFA
#   - subscale and total scores
#   - criterion correlations and binary gold-standard ROC analysis
#
# Optional package:
#   - lavaan, only when cfa = TRUE
#
# The returned object inherits from r4vn_tab, so the current tabexport()
# can export it without a separate parser.
# =============================================================================

.r4vn_ts_or <- function(x, y) if (is.null(x) || !length(x)) y else x

.r4vn_ts_stop <- function(...) stop(..., call. = FALSE)
.r4vn_ts_warn <- function(...) warning(..., call. = FALSE, immediate. = TRUE)

.r4vn_ts_escape <- function(x) {
  x <- as.character(x)
  x[is.na(x)] <- ""
  x <- gsub("&", "&amp;", x, fixed = TRUE)
  x <- gsub("<", "&lt;", x, fixed = TRUE)
  x <- gsub(">", "&gt;", x, fixed = TRUE)
  x <- gsub('"', "&quot;", x, fixed = TRUE)
  x
}

.r4vn_ts_num <- function(x, digit = 2L) {
  out <- rep("", length(x))
  ok <- is.finite(x)
  out[ok] <- formatC(x[ok], format = "f", digits = digit)
  out[is.infinite(x) & x > 0] <- "Inf"
  out[is.infinite(x) & x < 0] <- "-Inf"
  out
}

.r4vn_ts_safe <- function(x, fun, ...) {
  x <- x[is.finite(x)]
  if (!length(x)) return(NA_real_)
  fun(x, ...)
}

.r4vn_ts_p <- function(x, digit = 3L) {
  out <- rep("", length(x))
  ok <- is.finite(x)
  cut <- 10^(-digit)
  out[ok & x < cut] <- paste0("<", formatC(cut, format = "f", digits = digit))
  out[ok & x >= cut] <- formatC(x[ok & x >= cut], format = "f", digits = digit)
  out
}

.r4vn_ts_label <- function(x, fallback) {
  lab <- attr(x, "label", exact = TRUE)
  if (is.null(lab) || !length(lab) || is.na(lab[1L]) || !nzchar(as.character(lab[1L]))) fallback else as.character(lab[1L])
}

.r4vn_ts_get_data <- function(data = NULL) {
  if (!is.null(data)) {
    if (!is.data.frame(data)) data <- as.data.frame(data)
    return(data)
  }

  resolver <- get0(".r4vn_resolve_analysis_data", mode = "function", inherits = TRUE)
  if (!is.null(resolver)) {
    z <- try(resolver(NULL), silent = TRUE)
    if (is.data.frame(z)) return(z)
  }

  use <- get0("usedf", mode = "function", inherits = TRUE)
  if (!is.null(use)) {
    z <- try(use(quiet = TRUE), silent = TRUE)
    if (is.data.frame(z)) return(z)
  }

  for (op in c("R4VN.active_data", "r4vn.active_data")) {
    z <- getOption(op)
    if (is.data.frame(z)) return(z)
  }

  for (en in c(".r4vn_env", ".R4VN_env", ".r4vn_state")) {
    e <- get0(en, inherits = TRUE)
    if (is.environment(e)) {
      for (nm in c("active_data", "active", "data")) {
        if (exists(nm, envir = e, inherits = FALSE)) {
          z <- get(nm, envir = e, inherits = FALSE)
          if (is.data.frame(z)) return(z)
        }
      }
    }
  }

  .r4vn_ts_stop("No active data frame. Supply `data` or run `usedf(data)` first.")
}

.r4vn_ts_commit_data <- function(data, data_expr, data_supplied, env) {
  commit <- get0(".r4vn_commit_context", mode = "function", inherits = TRUE)
  if (!is.null(commit)) {
    target <- NULL
    if (!isTRUE(data_supplied)) {
      active_target <- get0(".r4vn_active_target", mode = "function", inherits = TRUE)
      if (!is.null(active_target)) target <- try(active_target(), silent = TRUE)
    } else if (is.symbol(data_expr)) {
      object_target <- get0(".r4vn_object_target", mode = "function", inherits = TRUE)
      if (!is.null(object_target)) target <- try(object_target(data_expr, env), silent = TRUE)
    }
    if (!is.null(target) && !inherits(target, "try-error")) {
      done <- try(commit(target, data), silent = TRUE)
      if (!inherits(done, "try-error")) return(invisible(TRUE))
    }
  }

  if (isTRUE(data_supplied) && is.symbol(data_expr)) {
    assign(as.character(data_expr), data, envir = env)
    return(invisible(TRUE))
  }

  set_active <- get0(".r4vn_set_active", mode = "function", inherits = TRUE)
  if (!is.null(set_active)) {
    done <- try(set_active(data, name = "tabscale", quiet = TRUE), silent = TRUE)
    if (!inherits(done, "try-error")) return(invisible(TRUE))
  }

  use <- get0("usedf", mode = "function", inherits = TRUE)
  if (!is.null(use)) {
    done <- try(use(data, quiet = TRUE), silent = TRUE)
    if (!inherits(done, "try-error")) return(invisible(TRUE))
  }
  invisible(FALSE)
}

.r4vn_ts_names <- function(x, arg = "vars", allow_null = FALSE) {
  if (is.null(x)) {
    if (allow_null) return(character())
    .r4vn_ts_stop(sprintf("`%s` is required.", arg))
  }
  if (inherits(x, "r4vn_vars")) return(unique(as.character(x$variable)))
  if (is.data.frame(x) && "variable" %in% names(x)) return(unique(as.character(x$variable)))
  if (inherits(x, "formula")) return(unique(all.vars(x)))
  if (is.character(x)) return(unique(x[nzchar(x)]))
  if (is.symbol(x)) return(as.character(x))
  .r4vn_ts_stop(sprintf("`%s` must be created by vars(), a character vector, or a one-sided formula.", arg))
}

.r4vn_ts_expr_name <- function(expr, value, data_names, arg) {
  if (identical(expr, quote(NULL))) return(NULL)
  if (is.symbol(expr)) {
    nm <- as.character(expr)
    if (nm %in% data_names) return(nm)
  }
  if (inherits(value, "r4vn_vars")) {
    z <- unique(as.character(value$variable))
    if (length(z) == 1L) return(z)
  }
  if (is.character(value) && length(value) == 1L && value %in% data_names) return(value)
  .r4vn_ts_stop(sprintf("`%s` must identify one variable in data.", arg))
}

.r4vn_ts_factor_map <- function(factor, all_items = NULL) {
  if (is.null(factor)) return(list())
  if (!is.list(factor) || inherits(factor, "r4vn_vars")) {
    .r4vn_ts_stop("`factor` must be a named list, for example list(Physical = vars(q1, q2), Mental = vars(q3, q4)).")
  }
  if (!length(factor)) return(list())
  nms <- names(factor)
  if (is.null(nms)) nms <- rep("", length(factor))
  blank <- !nzchar(nms)
  nms[blank] <- paste0("Factor", which(blank))
  nms <- make.unique(nms)
  if ("Total" %in% nms) .r4vn_ts_stop("`Total` is reserved for the full scale; use another factor name.")
  out <- stats::setNames(vector("list", length(factor)), nms)
  for (i in seq_along(factor)) {
    out[[i]] <- .r4vn_ts_names(factor[[i]], sprintf("factor[[%s]]", nms[i]))
    if (length(out[[i]]) < 2L) .r4vn_ts_stop(sprintf("Factor `%s` must contain at least two items.", nms[i]))
  }
  if (!is.null(all_items)) {
    bad <- setdiff(unique(unlist(out, use.names = FALSE)), all_items)
    if (length(bad)) .r4vn_ts_stop("Factor items not found in `vars`: ", paste(bad, collapse = ", "))
  }
  out
}

.r4vn_ts_numeric_item <- function(x, variable) {
  labels <- NULL
  original <- x

  if (is.logical(x)) {
    value <- as.numeric(x)
    labels <- c("0" = "No", "1" = "Yes")
  } else if (is.factor(x)) {
    lv <- levels(x)
    suppressWarnings(num_lv <- as.numeric(lv))
    if (all(is.finite(num_lv))) {
      value <- num_lv[as.integer(x)]
      labels <- stats::setNames(lv, as.character(num_lv))
    } else {
      value <- as.numeric(x)
      labels <- stats::setNames(lv, as.character(seq_along(lv)))
    }
  } else if (is.character(x)) {
    suppressWarnings(num <- as.numeric(x))
    if (all(is.na(x) | is.finite(num))) {
      value <- num
    } else {
      lv <- unique(x[!is.na(x)])
      value <- match(x, lv)
      labels <- stats::setNames(lv, as.character(seq_along(lv)))
      .r4vn_ts_warn(sprintf("Character item `%s` was scored in first-observed order; use factor levels to define the intended order.", variable))
    }
  } else if (is.numeric(x) || is.integer(x)) {
    value <- as.numeric(x)
    lab_attr <- attr(original, "labels", exact = TRUE)
    if (!is.null(lab_attr) && length(lab_attr)) {
      nm <- names(lab_attr)
      if (is.null(nm)) nm <- as.character(lab_attr)
      labels <- stats::setNames(nm, as.character(as.numeric(lab_attr)))
    }
  } else {
    .r4vn_ts_stop(sprintf("Item `%s` must be numeric, logical, factor, ordered factor, or character.", variable))
  }

  list(value = value, labels = labels)
}

.r4vn_ts_near_pd <- function(R, eps = 1e-8) {
  dn <- dimnames(R)
  R <- (R + t(R)) / 2
  diag(R) <- 1
  e <- eigen(R, symmetric = TRUE)
  if (min(e$values) < eps) {
    vals <- pmax(e$values, eps)
    R <- e$vectors %*% diag(vals, nrow = length(vals)) %*% t(e$vectors)
    d <- sqrt(diag(R))
    R <- sweep(sweep(R, 1L, d, "/"), 2L, d, "/")
    diag(R) <- 1
  }
  dimnames(R) <- dn
  R
}

.r4vn_ts_cor <- function(X, method = c("pearson", "spearman"), use = c("pairwise", "complete")) {
  method <- match.arg(method)
  use <- match.arg(use)
  use_r <- if (use == "complete") "complete.obs" else "pairwise.complete.obs"
  R <- suppressWarnings(stats::cor(X, use = use_r, method = method))
  if (any(!is.finite(R))) .r4vn_ts_stop("The item correlation matrix contains undefined values. Check items with no variation or excessive missingness.")
  .r4vn_ts_near_pd(R)
}

.r4vn_ts_cov <- function(X, use = c("pairwise", "complete")) {
  use <- match.arg(use)
  use_r <- if (use == "complete") "complete.obs" else "pairwise.complete.obs"
  S <- suppressWarnings(stats::cov(X, use = use_r))
  if (any(!is.finite(S))) .r4vn_ts_stop("The item covariance matrix contains undefined values.")
  S
}

.r4vn_ts_alpha_from_cov <- function(S) {
  k <- ncol(S)
  if (k < 2L) return(NA_real_)
  den <- sum(S)
  if (!is.finite(den) || den <= 0) return(NA_real_)
  k / (k - 1) * (1 - sum(diag(S)) / den)
}

.r4vn_ts_alpha <- function(X, use = c("pairwise", "complete"), cor_method = "pearson", bootstrap = 0L, conf = 0.95) {
  use <- match.arg(use)
  k <- ncol(X)
  S <- .r4vn_ts_cov(X, use)
  R <- .r4vn_ts_cor(X, cor_method, use)
  alpha <- .r4vn_ts_alpha_from_cov(S)
  alpha_std <- .r4vn_ts_alpha_from_cov(R)

  item_total <- vapply(seq_len(k), function(j) {
    rest <- rowSums(X[, -j, drop = FALSE], na.rm = TRUE)
    nrest <- rowSums(!is.na(X[, -j, drop = FALSE]))
    rest[nrest == 0L] <- NA_real_
    suppressWarnings(stats::cor(X[, j], rest, use = "pairwise.complete.obs", method = cor_method))
  }, numeric(1))

  alpha_deleted <- vapply(seq_len(k), function(j) {
    if (k <= 2L) return(NA_real_)
    .r4vn_ts_alpha_from_cov(.r4vn_ts_cov(X[, -j, drop = FALSE], use))
  }, numeric(1))

  off <- R[lower.tri(R)]
  mean_inter <- if (length(off)) mean(off, na.rm = TRUE) else NA_real_
  median_inter <- if (length(off)) stats::median(off, na.rm = TRUE) else NA_real_

  split_coefficient <- function(left) {
    right <- setdiff(seq_len(k), left)
    if (!length(left) || !length(right)) return(c(r = NA_real_, sb = NA_real_))
    a <- rowMeans(X[, left, drop = FALSE], na.rm = TRUE)
    b <- rowMeans(X[, right, drop = FALSE], na.rm = TRUE)
    a[rowSums(!is.na(X[, left, drop = FALSE])) == 0L] <- NA_real_
    b[rowSums(!is.na(X[, right, drop = FALSE])) == 0L] <- NA_real_
    r <- suppressWarnings(stats::cor(a, b, use = "pairwise.complete.obs", method = cor_method))
    sb <- if (is.finite(r) && r > -1) 2 * r / (1 + r) else NA_real_
    c(r = r, sb = sb)
  }
  odd <- seq(1L, k, by = 2L)
  split_odd_even <- split_coefficient(odd)
  split_values <- split_odd_even["sb"]
  if (k >= 4L) {
    size <- floor(k / 2L)
    if (choose(k, size) <= 500L) {
      candidates <- utils::combn(k, size, simplify = FALSE)
    } else {
      candidates <- replicate(500L, sort(sample.int(k, size)), simplify = FALSE)
    }
    split_values <- vapply(candidates, function(left) split_coefficient(left)["sb"], numeric(1))
  }
  split_values <- split_values[is.finite(split_values)]
  split_mean <- if (length(split_values)) mean(split_values) else NA_real_
  split_best <- if (length(split_values)) max(split_values) else NA_real_

  invR <- try(solve(R), silent = TRUE)
  lambda6 <- NA_real_
  if (!inherits(invR, "try-error")) {
    smc <- pmax(0, pmin(1, 1 - 1 / diag(invR)))
    error_var <- diag(S) * (1 - smc)
    total_var <- sum(S)
    if (is.finite(total_var) && total_var > 0) lambda6 <- 1 - sum(error_var) / total_var
  }

  total_var <- sum(S)
  lambda1 <- if (is.finite(total_var) && total_var > 0) 1 - sum(diag(S)) / total_var else NA_real_
  off_sq <- S^2
  diag(off_sq) <- 0
  lambda2 <- if (is.finite(lambda1) && total_var > 0) {
    lambda1 + sqrt(k / (k - 1) * sum(off_sq)) / total_var
  } else NA_real_
  lambda3 <- alpha
  lambda4 <- split_best
  row_off_sq <- rowSums(off_sq)
  lambda5 <- if (is.finite(lambda1) && total_var > 0) {
    lambda1 + 2 * sqrt(k / (k - 1)) * sqrt(max(row_off_sq)) / total_var
  } else NA_real_

  unique_values <- lapply(X, function(z) unique(z[is.finite(z)]))
  binary_items <- all(vapply(unique_values, length, integer(1)) <= 2L)
  kr20 <- if (binary_items) alpha else NA_real_

  omega <- NA_real_
  if (k >= 3L) {
    pa <- try(.r4vn_ts_paf(R, 1L, rotation = "none"), silent = TRUE)
    if (!inherits(pa, "try-error")) {
      lam <- as.numeric(pa$loadings[, 1L])
      if (sum(lam, na.rm = TRUE) < 0) lam <- -lam
      uniq <- pmax(0, 1 - lam^2)
      den <- sum(lam)^2 + sum(uniq)
      if (is.finite(den) && den > 0) omega <- sum(lam)^2 / den
    }
  }

  n_complete <- sum(stats::complete.cases(X))
  n_ci <- if (n_complete >= 4L) n_complete else nrow(X)
  ci <- c(lower = NA_real_, upper = NA_real_)
  omega_ci <- c(lower = NA_real_, upper = NA_real_)
  ci_method <- "Feldt"
  if (is.finite(alpha) && n_ci >= 4L && k >= 2L) {
    a <- 1 - conf
    ci <- c(
      lower = 1 - (1 - alpha) * stats::qf(1 - a / 2, n_ci - 1, (n_ci - 1) * (k - 1)),
      upper = 1 - (1 - alpha) * stats::qf(a / 2, n_ci - 1, (n_ci - 1) * (k - 1))
    )
    ci <- pmax(-1, pmin(1, ci))
  }
  bootstrap <- as.integer(bootstrap)
  if (bootstrap > 0L && nrow(X) >= 20L) {
    vals <- replicate(bootstrap, {
      id <- sample.int(nrow(X), nrow(X), replace = TRUE)
      tryCatch(.r4vn_ts_alpha_from_cov(.r4vn_ts_cov(X[id, , drop = FALSE], use)), error = function(e) NA_real_)
    })
    vals <- vals[is.finite(vals)]
    if (length(vals) >= max(20L, floor(bootstrap * 0.5))) {
      ci <- stats::quantile(vals, probs = c((1 - conf) / 2, 1 - (1 - conf) / 2), na.rm = TRUE, names = FALSE)
      names(ci) <- c("lower", "upper")
      ci_method <- "Bootstrap percentile"
    }
    omega_vals <- replicate(bootstrap, {
      id <- sample.int(nrow(X), nrow(X), replace = TRUE)
      tryCatch({
        rb <- .r4vn_ts_cor(X[id, , drop = FALSE], cor_method, use)
        pa <- .r4vn_ts_paf(rb, 1L, rotation = "none")
        lam <- as.numeric(pa$loadings[,1L]); if (sum(lam,na.rm=TRUE)<0) lam <- -lam
        uniq <- pmax(0,1-lam^2); sum(lam)^2/(sum(lam)^2+sum(uniq))
      }, error=function(e) NA_real_)
    })
    omega_vals <- omega_vals[is.finite(omega_vals)]
    if (length(omega_vals) >= max(20L, floor(bootstrap*.5))) {
      omega_ci <- stats::quantile(omega_vals, probs=c((1-conf)/2,1-(1-conf)/2),na.rm=TRUE,names=FALSE)
      names(omega_ci) <- c("lower","upper")
    }
  }

  ordinal_alpha <- omega_hierarchical <- glb <- NA_real_
  polychoric <- NULL
  if (requireNamespace("psych", quietly = TRUE) && k >= 3L &&
      all(vapply(X, function(z) length(unique(z[is.finite(z)])) <= 10L, logical(1)))) {
    pc <- try(psych::polychoric(X)$rho, silent = TRUE)
    if (!inherits(pc, "try-error") && is.matrix(pc) && all(is.finite(pc))) {
      polychoric <- .r4vn_ts_near_pd(pc)
      ordinal_alpha <- .r4vn_ts_alpha_from_cov(polychoric)
      om <- try(suppressMessages(psych::omega(polychoric, nfactors = min(3L, k - 1L),
                                               n.obs = n_ci, plot = FALSE)), silent = TRUE)
      if (!inherits(om, "try-error")) {
        if (length(om$omega_h) && is.finite(om$omega_h[1L])) omega_hierarchical <- om$omega_h[1L]
      }
    }
  }
  if (requireNamespace("psych", quietly = TRUE)) {
    glb_fit <- try(psych::glb.algebraic(S), silent = TRUE)
    if (!inherits(glb_fit, "try-error") && length(glb_fit$glb) && is.finite(glb_fit$glb[1L])) glb <- glb_fit$glb[1L]
  }

  list(
    alpha = alpha,
    alpha_std = alpha_std,
    alpha_ci = ci,
    alpha_ci_method = ci_method,
    ordinal_alpha = ordinal_alpha,
    item_total = item_total,
    alpha_deleted = alpha_deleted,
    mean_interitem = mean_inter,
    median_interitem = median_inter,
    split_half_r = unname(split_odd_even["r"]),
    spearman_brown = unname(split_odd_even["sb"]),
    split_half_mean = split_mean,
    split_half_best = split_best,
    lambda1 = lambda1,
    lambda2 = lambda2,
    lambda3 = lambda3,
    lambda4 = lambda4,
    lambda5 = lambda5,
    lambda6 = lambda6,
    glb = glb,
    kr20 = kr20,
    omega_total = omega,
    omega_ci = omega_ci,
    omega_hierarchical = omega_hierarchical,
    polychoric = polychoric,
    covariance = S,
    correlation = R
  )
}

.r4vn_ts_paf <- function(R, factors = 1L, rotation = c("varimax", "promax", "none"), maxit = 100L, tol = 1e-5) {
  rotation <- match.arg(rotation)
  p <- ncol(R)
  factors <- max(1L, min(as.integer(factors), p - 1L))
  invR <- try(solve(R), silent = TRUE)
  h2 <- if (inherits(invR, "try-error")) rep(0.5, p) else pmax(0.05, pmin(0.99, 1 - 1 / diag(invR)))
  L <- matrix(0, p, factors)

  for (i in seq_len(maxit)) {
    reduced <- R
    diag(reduced) <- h2
    e <- eigen(reduced, symmetric = TRUE)
    vals <- pmax(e$values[seq_len(factors)], 0)
    Lnew <- e$vectors[, seq_len(factors), drop = FALSE] %*% diag(sqrt(vals), nrow = factors)
    hnew <- rowSums(Lnew^2)
    if (max(abs(hnew - h2), na.rm = TRUE) < tol) {
      L <- Lnew
      h2 <- hnew
      break
    }
    L <- Lnew
    h2 <- pmin(0.999, pmax(0.001, hnew))
  }

  phi <- diag(factors)
  rotmat <- diag(factors)
  if (factors > 1L && rotation == "varimax") {
    vr <- stats::varimax(L)
    L <- unclass(vr$loadings)
    rotmat <- vr$rotmat
  } else if (factors > 1L && rotation == "promax") {
    pr <- stats::promax(L)
    L <- unclass(pr$loadings)
    rotmat <- pr$rotmat
    phi_try <- try(solve(t(rotmat) %*% rotmat), silent = TRUE)
    if (!inherits(phi_try, "try-error")) phi <- phi_try
  }

  colnames(L) <- paste0("F", seq_len(factors))
  rownames(L) <- rownames(R)
  dimnames(phi) <- list(colnames(L), colnames(L))
  list(loadings = L, communality = h2, uniqueness = pmax(0, 1 - h2), phi = phi, rotmat = rotmat, method = "Principal axis")
}

.r4vn_ts_kmo_bartlett <- function(R, n) {
  p <- ncol(R)
  invR <- try(solve(R), silent = TRUE)
  if (inherits(invR, "try-error")) {
    msa <- stats::setNames(rep(NA_real_, p), colnames(R))
    return(list(kmo = NA_real_, msa = msa, bartlett = c(chisq = NA, df = p * (p - 1) / 2, p = NA), determinant = NA_real_))
  }
  partial <- -stats::cov2cor(invR)
  diag(partial) <- 0
  r2 <- R^2
  diag(r2) <- 0
  p2 <- partial^2
  sum_r2 <- sum(r2) / 2
  sum_p2 <- sum(p2) / 2
  kmo <- sum_r2 / (sum_r2 + sum_p2)
  msa <- rowSums(r2) / (rowSums(r2) + rowSums(p2))
  names(msa) <- colnames(R)

  detR <- determinant(R, logarithm = TRUE)
  logdet <- as.numeric(detR$modulus)
  chi <- -(n - 1 - (2 * p + 5) / 6) * logdet
  df <- p * (p - 1) / 2
  pv <- stats::pchisq(chi, df = df, lower.tail = FALSE)
  list(kmo = kmo, msa = msa, bartlett = c(chisq = chi, df = df, p = pv), determinant = exp(logdet))
}

.r4vn_ts_parallel <- function(X, R, iter = 100L, quantile = 0.95, seed = NULL) {
  iter <- max(20L, as.integer(iter))
  n <- sum(stats::complete.cases(X))
  if (n < 20L) n <- nrow(X)
  p <- ncol(X)
  obs <- eigen(R, symmetric = TRUE, only.values = TRUE)$values
  if (!is.null(seed)) {
    old <- if (exists(".Random.seed", envir = .GlobalEnv, inherits = FALSE)) get(".Random.seed", envir = .GlobalEnv) else NULL
    on.exit({
      if (is.null(old)) {
        if (exists(".Random.seed", envir = .GlobalEnv, inherits = FALSE)) rm(".Random.seed", envir = .GlobalEnv)
      } else assign(".Random.seed", old, envir = .GlobalEnv)
    }, add = TRUE)
    set.seed(seed)
  }
  sim <- replicate(iter, {
    z <- matrix(stats::rnorm(n * p), nrow = n, ncol = p)
    eigen(stats::cor(z), symmetric = TRUE, only.values = TRUE)$values
  })
  random_q <- apply(sim, 1L, stats::quantile, probs = quantile, na.rm = TRUE)
  suggested <- sum(obs > random_q)
  if (suggested < 1L) suggested <- 1L
  tab <- data.frame(Component = seq_len(p), Observed = obs, Random95 = random_q, Retain = obs > random_q, stringsAsFactors = FALSE)
  list(table = tab, suggested = suggested)
}

.r4vn_ts_complexity <- function(L) {
  sq <- L^2
  den <- rowSums(sq^2)
  num <- rowSums(sq)^2
  ifelse(den > 0, num / den, NA_real_)
}

.r4vn_ts_efa <- function(X, cor_method, use, method, rotation, nfactor, parallel_iter, seed) {
  R <- .r4vn_ts_cor(X, cor_method, use)
  complete_n <- sum(stats::complete.cases(X))
  n_for_tests <- if (complete_n >= 20L) complete_n else nrow(X)
  diagnostics <- .r4vn_ts_kmo_bartlett(R, n_for_tests)
  parallel <- .r4vn_ts_parallel(X, R, parallel_iter, seed = seed)
  nf <- .r4vn_ts_or(nfactor, parallel$suggested)
  nf <- max(1L, min(as.integer(nf), ncol(X) - 1L))

  fit <- NULL
  if (method == "ml") {
    fa <- try(stats::factanal(covmat = R, factors = nf, n.obs = n_for_tests, rotation = rotation), silent = TRUE)
    if (!inherits(fa, "try-error")) {
      L <- unclass(fa$loadings)
      colnames(L) <- paste0("F", seq_len(ncol(L)))
      rownames(L) <- colnames(X)
      phi <- diag(ncol(L))
      if (rotation == "promax" && !is.null(fa$rotmat)) {
        phi_try <- try(solve(t(fa$rotmat) %*% fa$rotmat), silent = TRUE)
        if (!inherits(phi_try, "try-error")) phi <- phi_try
      }
      dimnames(phi) <- list(colnames(L), colnames(L))
      fit <- list(
        loadings = L,
        communality = 1 - fa$uniquenesses,
        uniqueness = fa$uniquenesses,
        phi = phi,
        method = "Maximum likelihood",
        statistic = unname(.r4vn_ts_or(fa$STATISTIC, NA_real_)),
        df = unname(.r4vn_ts_or(fa$dof, NA_real_)),
        p = unname(.r4vn_ts_or(fa$PVAL, NA_real_))
      )
    } else {
      .r4vn_ts_warn("Maximum-likelihood EFA did not converge; principal-axis factoring was used instead.")
    }
  }
  if (is.null(fit)) {
    fit <- .r4vn_ts_paf(R, nf, rotation)
    fit$statistic <- fit$df <- fit$p <- NA_real_
  }

  L <- fit$loadings
  dominant <- apply(abs(L), 1L, function(z) colnames(L)[which.max(z)])
  dominant_loading <- apply(L, 1L, function(z) z[which.max(abs(z))])
  loading_table <- data.frame(
    Variable = rownames(L),
    L,
    DominantFactor = dominant,
    DominantLoading = dominant_loading,
    Communality = fit$communality,
    Uniqueness = fit$uniqueness,
    Complexity = .r4vn_ts_complexity(L),
    check.names = FALSE,
    stringsAsFactors = FALSE
  )

  list(
    R = R,
    diagnostics = diagnostics,
    parallel = parallel,
    nfactor = nf,
    fit = fit,
    loadings = loading_table
  )
}

.r4vn_ts_score_one <- function(X, method = c("mean", "sum"), min_valid = NULL) {
  method <- match.arg(method)
  k <- ncol(X)
  if (is.null(min_valid)) min_valid <- ceiling(0.8 * k)
  if (length(min_valid) != 1L || !is.finite(min_valid)) .r4vn_ts_stop("`min_valid` must be one finite number.")
  if (min_valid > 0 && min_valid <= 1) min_valid <- ceiling(min_valid * k)
  min_valid <- max(1L, min(k, as.integer(min_valid)))
  nvalid <- rowSums(!is.na(X))
  score <- if (method == "mean") rowMeans(X, na.rm = TRUE) else rowSums(X, na.rm = TRUE)
  score[nvalid < min_valid] <- NA_real_
  list(score = score, nvalid = nvalid, min_valid = min_valid)
}

.r4vn_ts_subscale_scores <- function(X, factor_map, method, min_valid) {
  if (length(min_valid) > 1L && (is.null(names(min_valid)) || any(!nzchar(names(min_valid))))) {
    .r4vn_ts_stop("When `min_valid` has several values, it must be named by factor and/or `Total`.")
  }
  choose_min <- function(nm) {
    if (is.null(min_valid)) return(NULL)
    if (length(min_valid) == 1L) return(min_valid)
    if (nm %in% names(min_valid)) return(unname(min_valid[[nm]]))
    NULL
  }
  out <- list()
  if (length(factor_map)) {
    for (nm in names(factor_map)) out[[nm]] <- .r4vn_ts_score_one(X[, factor_map[[nm]], drop = FALSE], method, choose_min(nm))
  }
  out[["Total"]] <- .r4vn_ts_score_one(X, method, choose_min("Total"))
  out
}

.r4vn_ts_roc <- function(y, score, event = NULL, direction = c("auto", "higher", "lower"), conf = 0.95) {
  direction <- match.arg(direction)
  ok <- !is.na(y) & is.finite(score)
  y <- y[ok]
  score <- score[ok]
  if (is.factor(y)) y <- droplevels(y)
  if (is.factor(y)) {
    lv <- levels(y)
    if (length(lv) != 2L) .r4vn_ts_stop("`gold` must have exactly two levels for ROC analysis.")
    event <- .r4vn_ts_or(event, lv[2L])
    if (!event %in% lv) .r4vn_ts_stop("`event` was not found among gold-standard levels.")
    yy <- as.integer(y == event)
    nonevent <- lv[lv != event][1L]
  } else {
    vals <- sort(unique(y[!is.na(y)]))
    if (length(vals) != 2L) .r4vn_ts_stop("`gold` must have exactly two values for ROC analysis.")
    event <- .r4vn_ts_or(event, vals[2L])
    yy <- as.integer(y == event)
    nonevent <- vals[vals != event][1L]
  }
  n1 <- sum(yy == 1L)
  n0 <- sum(yy == 0L)
  if (n1 < 2L || n0 < 2L) .r4vn_ts_stop("ROC analysis requires at least two event and two non-event observations.")

  auc_higher <- (sum(rank(score, ties.method = "average")[yy == 1L]) - n1 * (n1 + 1) / 2) / (n1 * n0)
  if (direction == "auto") direction <- if (auc_higher >= 0.5) "higher" else "lower"
  auc <- if (direction == "higher") auc_higher else 1 - auc_higher

  u <- sort(unique(score))
  cuts <- if (length(u) == 1L) u else c(u[1L] - .Machine$double.eps, (u[-1L] + u[-length(u)]) / 2, u[length(u)] + .Machine$double.eps)
  rows <- lapply(cuts, function(cut) {
    pred <- if (direction == "higher") score >= cut else score <= cut
    tp <- sum(pred & yy == 1L)
    fn <- sum(!pred & yy == 1L)
    tn <- sum(!pred & yy == 0L)
    fp <- sum(pred & yy == 0L)
    sens <- tp / (tp + fn)
    spec <- tn / (tn + fp)
    ppv <- if (tp + fp > 0) tp / (tp + fp) else NA_real_
    npv <- if (tn + fn > 0) tn / (tn + fn) else NA_real_
    lrpos <- if (1 - spec > 0) sens / (1 - spec) else Inf
    lrneg <- if (spec > 0) (1 - sens) / spec else Inf
    accuracy <- (tp + tn) / (tp + tn + fp + fn)
    dor <- if (is.finite(lrpos) && is.finite(lrneg) && lrneg > 0) lrpos / lrneg else Inf
    c(cutoff = cut, sensitivity = sens, specificity = spec, ppv = ppv, npv = npv, lr_pos = lrpos, lr_neg = lrneg, accuracy = accuracy, dor = dor, youden = sens + spec - 1, tp = tp, fp = fp, tn = tn, fn = fn)
  })
  curve <- as.data.frame(do.call(rbind, rows), stringsAsFactors = FALSE)
  best <- curve[which.max(curve$youden), , drop = FALSE]

  q1 <- auc / (2 - auc)
  q2 <- 2 * auc^2 / (1 + auc)
  se <- sqrt((auc * (1 - auc) + (n1 - 1) * (q1 - auc^2) + (n0 - 1) * (q2 - auc^2)) / (n1 * n0))
  z <- stats::qnorm(1 - (1 - conf) / 2)
  ci <- pmax(0, pmin(1, auc + c(-1, 1) * z * se))
  auc_z <- if (is.finite(se) && se > 0) (auc - 0.5) / se else NA_real_
  auc_p <- if (is.finite(auc_z)) 2 * stats::pnorm(abs(auc_z), lower.tail = FALSE) else NA_real_

  list(
    event = as.character(event),
    nonevent = as.character(nonevent),
    direction = direction,
    n_event = n1,
    n_nonevent = n0,
    auc = auc,
    se = se,
    p = auc_p,
    ci = stats::setNames(ci, c("lower", "upper")),
    best = best,
    curve = curve[order(1 - curve$specificity, curve$sensitivity), , drop = FALSE]
  )
}

.r4vn_ts_htmt <- function(R, factor_map) {
  if (length(factor_map) < 2L) return(data.frame())
  out <- list()
  nms <- names(factor_map)
  z <- 0L
  for (i in seq_len(length(nms) - 1L)) {
    for (j in (i + 1L):length(nms)) {
      a <- factor_map[[i]]
      b <- factor_map[[j]]
      cross <- mean(abs(R[a, b, drop = FALSE]), na.rm = TRUE)
      wa <- if (length(a) > 1L) mean(abs(R[a, a, drop = FALSE][lower.tri(R[a, a, drop = FALSE])]), na.rm = TRUE) else NA_real_
      wb <- if (length(b) > 1L) mean(abs(R[b, b, drop = FALSE][lower.tri(R[b, b, drop = FALSE])]), na.rm = TRUE) else NA_real_
      htmt <- cross / sqrt(wa * wb)
      z <- z + 1L
      out[[z]] <- data.frame(Factor1 = nms[i], Factor2 = nms[j], HTMT = htmt, stringsAsFactors = FALSE)
    }
  }
  do.call(rbind, out)
}

.r4vn_ts_cfa <- function(X, original_data, factor_map, ordered, estimator, missing, cfa_group_name, modification, modification_min) {
  if (!requireNamespace("lavaan", quietly = TRUE)) {
    .r4vn_ts_stop("CFA requires the optional package `lavaan`. Install it with install.packages(\"lavaan\") and add lavaan to Suggests in DESCRIPTION.")
  }
  if (!length(factor_map)) .r4vn_ts_stop("`factor` must be declared when `cfa = TRUE`.")

  item_labels <- colnames(X)
  item_model_names <- make.names(item_labels, unique = TRUE)
  item_to_model <- stats::setNames(item_model_names, item_labels)
  model_to_item <- stats::setNames(item_labels, item_model_names)

  factor_labels <- names(factor_map)
  factor_model_names <- make.names(factor_labels, unique = TRUE)
  factor_to_model <- stats::setNames(factor_model_names, factor_labels)
  model_to_factor <- stats::setNames(factor_labels, factor_model_names)

  factor_map_model <- lapply(factor_map, function(items) unname(item_to_model[items]))
  names(factor_map_model) <- factor_model_names
  model <- paste(vapply(names(factor_map_model), function(nm) paste0(nm, " =~ ", paste(factor_map_model[[nm]], collapse = " + ")), character(1)), collapse = "\n")
  item_display <- stats::setNames(vapply(item_labels, function(v) .r4vn_ts_label(original_data[[v]], v), character(1)), item_labels)
  display_model <- paste(vapply(factor_labels, function(nm) paste0(nm, " =~ ", paste(item_display[factor_map[[nm]]], collapse = " + ")), character(1)), collapse = "\n")

  cfa_data <- as.data.frame(X, check.names = FALSE)
  names(cfa_data) <- item_model_names
  group_model_name <- NULL
  if (!is.null(cfa_group_name)) {
    group_model_name <- make.unique(c(names(cfa_data), ".r4vn_cfa_group"), sep = "_")
    group_model_name <- tail(group_model_name, 1L)
    cfa_data[[group_model_name]] <- original_data[[cfa_group_name]]
  }

  ordered_labels <- character()
  if (isTRUE(ordered)) ordered_labels <- item_labels
  else if (inherits(ordered, "r4vn_vars")) ordered_labels <- intersect(as.character(ordered$variable), item_labels)
  else if (is.character(ordered)) ordered_labels <- intersect(ordered, item_labels)
  ordered_names <- unname(item_to_model[ordered_labels])

  estimator_use <- estimator
  if (identical(estimator, "auto")) estimator_use <- if (length(ordered_names)) "WLSMV" else "MLR"

  args <- list(model = model, data = cfa_data, estimator = estimator_use, std.lv = TRUE)
  if (length(ordered_names)) args$ordered <- ordered_names
  if (!is.null(group_model_name)) args$group <- group_model_name
  if (!length(ordered_names) && !is.null(missing) && nzchar(missing)) args$missing <- missing

  fit <- try(do.call(lavaan::cfa, args), silent = TRUE)
  if (inherits(fit, "try-error")) .r4vn_ts_stop("CFA failed: ", as.character(fit))

  measures_wanted <- c("chisq", "df", "pvalue", "cfi", "tli", "rmsea", "rmsea.ci.lower", "rmsea.ci.upper", "srmr", "aic", "bic")
  measures <- try(lavaan::fitMeasures(fit, measures_wanted), silent = TRUE)
  if (inherits(measures, "try-error")) measures <- stats::setNames(rep(NA_real_, length(measures_wanted)), measures_wanted)

  std <- lavaan::standardizedSolution(fit)
  loading_index <- std$op == "=~"
  loading_columns <- c("lhs", "rhs", "est.std", "se", "z", "pvalue", "ci.lower", "ci.upper")
  if ("group" %in% names(std)) loading_columns <- c(loading_columns, "group")
  load <- std[loading_index, loading_columns, drop = FALSE]
  names(load)[seq_len(8L)] <- c("Factor", "Item", "Loading", "SE", "z", "p", "CI_lower", "CI_upper")

  group_labels <- NULL
  if ("group" %in% names(load)) {
    group_labels <- try(lavaan::lavInspect(fit, "group.label"), silent = TRUE)
    if (inherits(group_labels, "try-error") || !length(group_labels)) group_labels <- sort(unique(load$group))
    load$Group <- as.character(group_labels[load$group])
    load$group <- NULL
  }

  r2 <- try(lavaan::inspect(fit, "r2"), silent = TRUE)
  if (!inherits(r2, "try-error")) {
    if (is.list(r2) && "Group" %in% names(load)) {
      group_index <- match(load$Group, as.character(group_labels))
      load$R2 <- vapply(seq_len(nrow(load)), function(i) unname(r2[[group_index[i]]][load$Item[i]]), numeric(1))
    } else {
      if (is.list(r2)) r2 <- r2[[1L]]
      load$R2 <- unname(r2[load$Item])
    }
  } else load$R2 <- NA_real_
  load$Factor <- unname(model_to_factor[load$Factor])
  load$Item <- unname(model_to_item[load$Item])
  load$`Item label` <- vapply(load$Item, function(v) .r4vn_ts_label(original_data[[v]], v), character(1))
  load <- load[, c(if ("Group" %in% names(load)) "Group", "Factor", "Item label", "Item",
                   setdiff(names(load), c("Group", "Factor", "Item label", "Item"))), drop = FALSE]
  if ("Group" %in% names(load)) load <- load[, c("Group", setdiff(names(load), "Group")), drop = FALSE]

  reliability_groups <- if ("Group" %in% names(load)) unique(load$Group) else NA_character_
  reliability <- do.call(rbind, lapply(reliability_groups, function(gr) {
    zgroup <- if (is.na(gr)) load else load[load$Group == gr, , drop = FALSE]
    do.call(rbind, lapply(factor_labels, function(nm) {
      z <- zgroup[zgroup$Factor == nm, , drop = FALSE]
      lam <- z$Loading
      theta <- pmax(0, 1 - lam^2)
      cr <- sum(lam)^2 / (sum(lam)^2 + sum(theta))
      ave <- sum(lam^2) / (sum(lam^2) + sum(theta))
      out <- data.frame(Factor = nm, CompositeReliability = cr, AVE = ave, SqrtAVE = sqrt(ave), stringsAsFactors = FALSE)
      if (!is.na(gr)) out <- data.frame(Group = gr, out, check.names = FALSE, stringsAsFactors = FALSE)
      out
    }))
  }))

  factor_cor_columns <- c("lhs", "rhs", "est.std", "pvalue")
  if ("group" %in% names(std)) factor_cor_columns <- c(factor_cor_columns, "group")
  factor_cor <- std[std$op == "~~" & std$lhs != std$rhs & std$lhs %in% factor_model_names & std$rhs %in% factor_model_names, factor_cor_columns, drop = FALSE]
  if (nrow(factor_cor)) {
    names(factor_cor)[seq_len(4L)] <- c("Factor1", "Factor2", "Correlation", "p")
    factor_cor$Factor1 <- unname(model_to_factor[factor_cor$Factor1])
    factor_cor$Factor2 <- unname(model_to_factor[factor_cor$Factor2])
    if ("group" %in% names(factor_cor)) {
      factor_cor$Group <- as.character(group_labels[factor_cor$group])
      factor_cor$group <- NULL
      factor_cor <- factor_cor[, c("Group", setdiff(names(factor_cor), "Group")), drop = FALSE]
    }
  }

  R <- .r4vn_ts_cor(X, "pearson", "pairwise")
  htmt <- .r4vn_ts_htmt(R, factor_map)

  mi <- data.frame()
  if (isTRUE(modification)) {
    mi0 <- try(lavaan::modindices(fit, sort. = TRUE), silent = TRUE)
    if (!inherits(mi0, "try-error")) {
      mi0 <- mi0[is.finite(mi0$mi) & mi0$mi >= modification_min, c("lhs", "op", "rhs", "mi", "epc", "sepc.all"), drop = FALSE]
      if (nrow(mi0) > 30L) mi0 <- mi0[seq_len(30L), , drop = FALSE]
      names(mi0) <- c("Left", "Operator", "Right", "MI", "EPC", "StandardizedEPC")
      mi0$Left <- ifelse(mi0$Left %in% names(model_to_factor), unname(model_to_factor[mi0$Left]),
                         ifelse(mi0$Left %in% names(model_to_item), unname(model_to_item[mi0$Left]), mi0$Left))
      mi0$Right <- ifelse(mi0$Right %in% names(model_to_factor), unname(model_to_factor[mi0$Right]),
                          ifelse(mi0$Right %in% names(model_to_item), unname(model_to_item[mi0$Right]), mi0$Right))
      mi <- mi0
    }
  }

  list(
    fit = fit,
    model = display_model,
    model_syntax = model,
    estimator = estimator_use,
    ordered = ordered_labels,
    group = cfa_group_name,
    measures = measures,
    loadings = load,
    reliability = reliability,
    factor_correlations = factor_cor,
    htmt = htmt,
    modification = mi
  )
}

.r4vn_ts_bind_fill <- function(tables) {
  tables <- tables[vapply(tables, function(x) is.data.frame(x) && nrow(x) > 0L, logical(1))]
  if (!length(tables)) return(data.frame())
  all_names <- unique(unlist(lapply(tables, names), use.names = FALSE))
  out <- lapply(names(tables), function(section) {
    x <- tables[[section]]
    miss <- setdiff(all_names, names(x))
    for (nm in miss) x[[nm]] <- ""
    x <- x[, all_names, drop = FALSE]
    data.frame(Section = section, x, check.names = FALSE, stringsAsFactors = FALSE)
  })
  do.call(rbind, out)
}

.r4vn_ts_transpose_display <- function(x) {
  if (!is.data.frame(x) || !nrow(x) || !ncol(x)) return(x)

  # A single result is clearest as Statistic | Value.
  if (nrow(x) == 1L) {
    return(data.frame(
      Statistic = names(x),
      Value = unlist(x[1L, , drop = FALSE], use.names = FALSE),
      check.names = FALSE, stringsAsFactors = FALSE
    ))
  }

  # Prefer a meaningful row identifier so scales/models become columns.
  identifiers <- c("Scale", "Model", "Factor", "Group", "Outcome", "Item", "Variable")
  id <- identifiers[identifiers %in% names(x)]
  id <- id[vapply(id, function(nm) {
    z <- as.character(x[[nm]])
    all(!is.na(z) & nzchar(z)) && !anyDuplicated(z)
  }, logical(1))]
  id <- if (length(id)) id[1L] else NULL

  if (is.null(id)) {
    labels <- paste("Result", seq_len(nrow(x)))
    value_names <- names(x)
  } else {
    labels <- make.unique(as.character(x[[id]]))
    # `Analysis` repeats the section caption and need not consume a table row.
    value_names <- setdiff(names(x), c(id, "Analysis"))
  }

  out <- data.frame(Statistic = value_names, check.names = FALSE,
                    stringsAsFactors = FALSE)
  for (i in seq_len(nrow(x))) {
    out[[labels[i]]] <- unlist(x[i, value_names, drop = FALSE], use.names = FALSE)
  }
  out
}

.r4vn_ts_html_table <- function(df, caption = NULL, digits = 3L,
                                 p_digits = digits, transpose = FALSE) {
  if (!is.data.frame(df) || !nrow(df)) return("")
  x <- df
  for (j in seq_along(x)) {
    if (is.numeric(x[[j]])) {
      is_p <- grepl("(^p$|(^|[ _.])p($|[ _.])|pvalue|p-value|p value)", tolower(names(x)[j]))
      x[[j]] <- if (is_p) .r4vn_ts_p(x[[j]], p_digits) else .r4vn_ts_num(x[[j]], digits)
    } else x[[j]] <- as.character(x[[j]])
  }
  if (isTRUE(transpose)) x <- .r4vn_ts_transpose_display(x)
  head <- paste0("<tr>", paste0("<th>", .r4vn_ts_escape(names(x)), "</th>", collapse = ""), "</tr>")
  body <- vapply(seq_len(nrow(x)), function(i) {
    paste0("<tr>", paste0("<td>", .r4vn_ts_escape(unlist(x[i, , drop = FALSE], use.names = FALSE)), "</td>", collapse = ""), "</tr>")
  }, character(1))
  cap <- if (!is.null(caption) && nzchar(caption)) paste0("<h3>", .r4vn_ts_escape(caption), "</h3>") else ""
  paste0("<div class=\"ts-section\">", cap, "<table><thead>", head, "</thead><tbody>", paste(body, collapse = ""), "</tbody></table></div>")
}

.r4vn_ts_svg_line <- function(x, ys, labels, width = 720, height = 330, xlab = "", ylab = "") {
  if (!length(x) || !length(ys)) return("")
  ys <- lapply(ys, as.numeric)
  xr <- range(x, finite = TRUE)
  yr <- range(unlist(ys), finite = TRUE)
  if (!all(is.finite(c(xr, yr)))) return("")
  if (diff(xr) == 0) xr <- xr + c(-0.5, 0.5)
  if (diff(yr) == 0) yr <- yr + c(-0.5, 0.5)
  px <- function(z) 55 + (z - xr[1]) / diff(xr) * (width - 80)
  py <- function(z) height - 40 - (z - yr[1]) / diff(yr) * (height - 70)
  paths <- vapply(seq_along(ys), function(i) {
    y <- ys[[i]]
    ok <- is.finite(x) & is.finite(y)
    pts <- paste0(round(px(x[ok]), 1), ",", round(py(y[ok]), 1), collapse = " ")
    dash <- if (i == 1L) "" else " stroke-dasharray=\"6,4\""
    paste0("<polyline fill=\"none\" stroke=\"currentColor\" stroke-width=\"2\"", dash, " points=\"", pts, "\"/>")
  }, character(1))
  legend <- paste0("<text x=\"", 60 + (seq_along(labels) - 1) * 160, "\" y=\"20\" font-size=\"12\">", .r4vn_ts_escape(labels), "</text>", collapse = "")
  paste0(
    "<svg class=\"ts-svg\" viewBox=\"0 0 ", width, " ", height, "\" role=\"img\">",
    "<line x1=\"55\" y1=\"", height - 40, "\" x2=\"", width - 25, "\" y2=\"", height - 40, "\" stroke=\"currentColor\"/>",
    "<line x1=\"55\" y1=\"30\" x2=\"55\" y2=\"", height - 40, "\" stroke=\"currentColor\"/>",
    paste(paths, collapse = ""), legend,
    "<text x=\"", width / 2, "\" y=\"", height - 8, "\" text-anchor=\"middle\" font-size=\"12\">", .r4vn_ts_escape(xlab), "</text>",
    "<text x=\"15\" y=\"", height / 2, "\" transform=\"rotate(-90 15 ", height / 2, ")\" text-anchor=\"middle\" font-size=\"12\">", .r4vn_ts_escape(ylab), "</text>",
    "</svg>"
  )
}

.r4vn_ts_clean_factor_name <- function(x) {
  z <- gsub("[^[:alnum:]_]+", "_", x)
  z <- gsub("^_+|_+$", "", z)
  ifelse(nzchar(z), z, "subscale")
}

.r4vn_ts_fisher_cor <- function(x, y, method = "pearson", conf = 0.95) {
  ok <- is.finite(x) & is.finite(y)
  n <- sum(ok)
  if (n < 4L) return(c(N = n, r = NA_real_, lower = NA_real_, upper = NA_real_, p = NA_real_))
  test <- try(suppressWarnings(stats::cor.test(x[ok], y[ok], method = method,
                                               exact = FALSE)), silent = TRUE)
  if (inherits(test, "try-error")) return(c(N = n, r = NA_real_, lower = NA_real_, upper = NA_real_, p = NA_real_))
  r <- unname(test$estimate)
  ci <- c(NA_real_, NA_real_)
  if (is.finite(r) && abs(r) < 1 && n > 3L) {
    z <- atanh(r)
    q <- stats::qnorm(1 - (1 - conf) / 2)
    ci <- tanh(z + c(-1, 1) * q / sqrt(n - 3))
  } else if (is.finite(r)) ci <- rep(r, 2L)
  c(N = n, r = r, lower = ci[1L], upper = ci[2L], p = test$p.value)
}

.r4vn_ts_external_validity <- function(scores, data, variables, hypothesis,
                                        method = "pearson", conf = 0.95,
                                        type = "Convergent") {
  if (!length(variables)) return(data.frame())
  bad <- setdiff(variables, names(data))
  if (length(bad)) .r4vn_ts_stop(type, " validity variables not found in data: ", paste(bad, collapse = ", "))
  out <- list()
  z <- 0L
  for (sn in names(scores)) {
    for (v in variables) {
      y <- .r4vn_ts_numeric_item(data[[v]], v)$value
      est <- .r4vn_ts_fisher_cor(scores[[sn]], y, method, conf)
      z <- z + 1L
      met <- if (is.null(hypothesis)) NA else {
        h <- hypothesis
        if (length(h) > 1L && !is.null(names(h)) && v %in% names(h)) h <- h[[v]]
        h <- as.numeric(h[1L])
        if (!is.finite(h)) NA else if (type == "Convergent") abs(est["r"]) >= abs(h) else abs(est["r"]) <= abs(h)
      }
      out[[z]] <- data.frame(
        Scale = sn,
        Variable = .r4vn_ts_label(data[[v]], v),
        `Variable name` = v,
        Type = type,
        N = unname(est["N"]),
        Method = tools::toTitleCase(method),
        Correlation = unname(est["r"]),
        `CI lower` = unname(est["lower"]),
        `CI upper` = unname(est["upper"]),
        p = unname(est["p"]),
        `Hypothesis threshold` = if (is.null(hypothesis)) NA_real_ else suppressWarnings(as.numeric(hypothesis[1L])),
        `Hypothesis met` = if (is.na(met)) "Not specified" else if (met) "Yes" else "No",
        check.names = FALSE, stringsAsFactors = FALSE
      )
    }
  }
  do.call(rbind, out)
}

.r4vn_ts_icc <- function(M, conf = 0.95) {
  M <- as.matrix(M)
  storage.mode(M) <- "double"
  M <- M[stats::complete.cases(M), , drop = FALSE]
  n <- nrow(M); k <- ncol(M)
  empty <- c(N = n, Raters = k, ICC_A1 = NA, ICC_Ak = NA, ICC_C1 = NA, ICC_Ck = NA,
             lower = NA, upper = NA, SEM = NA, MDC95 = NA)
  if (n < 3L || k < 2L) return(empty)
  grand <- mean(M)
  row_mean <- rowMeans(M)
  col_mean <- colMeans(M)
  ss_row <- k * sum((row_mean - grand)^2)
  ss_col <- n * sum((col_mean - grand)^2)
  ss_error <- sum((M - row_mean - rep(col_mean, each = n) + grand)^2)
  msr <- ss_row / (n - 1)
  msc <- ss_col / (k - 1)
  mse <- ss_error / ((n - 1) * (k - 1))
  icc_a1 <- (msr - mse) / (msr + (k - 1) * mse + k * (msc - mse) / n)
  icc_ak <- (msr - mse) / (msr + (msc - mse) / n)
  icc_c1 <- (msr - mse) / (msr + (k - 1) * mse)
  icc_ck <- (msr - mse) / msr
  f <- if (mse > 0) msr / mse else Inf
  a <- 1 - conf
  lower <- if (is.finite(f) && f > 0) (f / stats::qf(1 - a / 2, n - 1, (n - 1) * (k - 1)) - 1) /
    (f / stats::qf(1 - a / 2, n - 1, (n - 1) * (k - 1)) + k - 1) else NA_real_
  upper <- if (is.finite(f) && f > 0) (f / stats::qf(a / 2, n - 1, (n - 1) * (k - 1)) - 1) /
    (f / stats::qf(a / 2, n - 1, (n - 1) * (k - 1)) + k - 1) else NA_real_
  sem <- sqrt(max(mse, 0))
  c(N = n, Raters = k, ICC_A1 = icc_a1, ICC_Ak = icc_ak,
    ICC_C1 = icc_c1, ICC_Ck = icc_ck, lower = lower, upper = upper,
    SEM = sem, MDC95 = 1.96 * sqrt(2) * sem)
}

.r4vn_ts_kappa <- function(M, weighted = TRUE, conf = 0.95) {
  M <- as.data.frame(M, check.names = FALSE)
  M <- M[stats::complete.cases(M), , drop = FALSE]
  n <- nrow(M); k <- ncol(M)
  if (n < 3L || k < 2L) return(c(N = n, Raters = k, Agreement = NA, Kappa = NA, lower = NA, upper = NA))
  lev <- sort(unique(unlist(lapply(M, as.character), use.names = FALSE)))
  Z <- sapply(M, function(x) match(as.character(x), lev))
  if (k == 2L) {
    tab <- table(factor(Z[, 1], levels = seq_along(lev)), factor(Z[, 2], levels = seq_along(lev)))
    p <- tab / sum(tab)
    if (weighted && length(lev) > 2L) {
      w <- 1 - outer(seq_along(lev), seq_along(lev), function(a, b) (a - b)^2 / (length(lev) - 1)^2)
    } else w <- diag(length(lev))
    po <- sum(w * p)
    pe <- sum(w * outer(rowSums(p), colSums(p)))
    kap <- if (pe < 1) (po - pe) / (1 - pe) else NA_real_
    agreement <- sum(diag(tab)) / sum(tab)
  } else {
    counts <- t(apply(Z, 1L, function(z) tabulate(z, nbins = length(lev))))
    pi <- colSums(counts) / (n * k)
    pj <- (rowSums(counts^2) - k) / (k * (k - 1))
    po <- mean(pj); pe <- sum(pi^2)
    kap <- if (pe < 1) (po - pe) / (1 - pe) else NA_real_
    agreement <- mean(apply(counts, 1L, max) == k)
  }
  se <- if (is.finite(kap)) sqrt(max(po * (1 - po), 0) / (n * max((1 - pe)^2, .Machine$double.eps))) else NA_real_
  q <- stats::qnorm(1 - (1 - conf) / 2)
  ci <- if (is.finite(se)) pmax(-1, pmin(1, kap + c(-1, 1) * q * se)) else c(NA, NA)
  c(N = n, Raters = k, Agreement = agreement, Kappa = kap, lower = ci[1L], upper = ci[2L])
}

.r4vn_ts_agreement <- function(first, second, scale_name, type, conf = 0.95) {
  ok <- is.finite(first) & is.finite(second)
  x <- first[ok]; y <- second[ok]
  icc <- .r4vn_ts_icc(cbind(x, y), conf)
  cor_p <- .r4vn_ts_fisher_cor(x, y, "pearson", conf)
  cor_s <- .r4vn_ts_fisher_cor(x, y, "spearman", conf)
  d <- y - x
  md <- if (length(d)) mean(d) else NA_real_
  sdd <- if (length(d) > 1L) stats::sd(d) else NA_real_
  se_md <- sdd / sqrt(length(d))
  q <- if (length(d) > 1L) stats::qt(1 - (1 - conf) / 2, length(d) - 1L) else NA_real_
  tt <- if (length(d) > 1L) try(stats::t.test(y, x, paired = TRUE, conf.level = conf), silent = TRUE) else NULL
  wp <- if (length(d) > 1L) try(suppressWarnings(stats::wilcox.test(y, x, paired = TRUE, exact = FALSE)), silent = TRUE) else NULL
  pooled <- sqrt((stats::var(x) + stats::var(y)) / 2)
  grand_mean <- mean(c(x,y))
  cv_within <- if (is.finite(grand_mean) && abs(grand_mean) > .Machine$double.eps)
    100*sqrt(mean(d^2)/2)/abs(grand_mean) else NA_real_
  tab <- data.frame(
    Analysis = type, Scale = scale_name, N = length(d),
    `Time/form 1 mean` = mean(x), `Time/form 2 mean` = mean(y),
    `Mean difference` = md, `Difference CI lower` = md - q * se_md,
    `Difference CI upper` = md + q * se_md,
    `Paired t p` = if (!is.null(tt) && !inherits(tt, "try-error")) tt$p.value else NA_real_,
    `Wilcoxon p` = if (!is.null(wp) && !inherits(wp, "try-error")) wp$p.value else NA_real_,
    `Pearson r` = unname(cor_p["r"]), `Spearman rho` = unname(cor_s["r"]),
    `ICC absolute single` = unname(icc["ICC_A1"]),
    `ICC CI lower` = unname(icc["lower"]), `ICC CI upper` = unname(icc["upper"]),
    SEM = unname(icc["SEM"]), MDC95 = unname(icc["MDC95"]),
    `Within-subject CV %` = cv_within,
    `Cohen dz` = if (is.finite(sdd) && sdd > 0) md / sdd else NA_real_,
    `Standardized mean difference` = if (is.finite(pooled) && pooled > 0) md / pooled else NA_real_,
    `Bland-Altman lower` = md - 1.96 * sdd, `Bland-Altman upper` = md + 1.96 * sdd,
    check.names = FALSE, stringsAsFactors = FALSE
  )
  list(table = tab, first = x, second = y, difference = d)
}

.r4vn_ts_content_validity <- function(content, cutoff = 3, max_score = 4, conf = 0.95) {
  if (is.null(content)) return(NULL)
  if (!is.data.frame(content) && !is.matrix(content)) .r4vn_ts_stop("`content` must be an expert-by-item data frame or numeric matrix.")
  labs <- if (is.data.frame(content)) vapply(seq_along(content), function(j) .r4vn_ts_label(content[[j]], names(content)[j]), character(1)) else colnames(content)
  M <- as.matrix(content)
  suppressWarnings(storage.mode(M) <- "double")
  if (!ncol(M) || nrow(M) < 2L || !any(is.finite(M)) || any(!is.na(M) & !is.finite(M))) .r4vn_ts_stop("`content` must contain numeric expert relevance ratings.")
  if (!is.finite(cutoff) || !is.finite(max_score) || cutoff > max_score) .r4vn_ts_stop("`content_cutoff` and `content_max` must be finite, with cutoff not exceeding the maximum.")
  if (any(M > max_score, na.rm=TRUE)) .r4vn_ts_stop("Content-validity ratings exceed `content_max`.")
  if (is.null(colnames(M))) colnames(M) <- paste0("Item", seq_len(ncol(M)))
  if (is.null(labs) || length(labs) != ncol(M)) labs <- colnames(M)
  rows <- lapply(seq_len(ncol(M)), function(j) {
    z <- M[, j]; z <- z[is.finite(z)]
    n <- length(z); agree <- sum(z >= cutoff)
    icvi <- if (n) agree / n else NA_real_
    ci <- if (n) stats::binom.test(agree,n,conf.level=conf)$conf.int else c(NA_real_,NA_real_)
    pc <- if (n) choose(n, agree) * 0.5^n else NA_real_
    kstar <- if (is.finite(icvi) && is.finite(pc) && pc < 1) (icvi - pc) / (1 - pc) else NA_real_
    data.frame(Item = labs[j], Variable = colnames(M)[j], Experts = n,
               `Relevant ratings` = agree, `I-CVI` = icvi,
               `I-CVI CI lower` = ci[1L], `I-CVI CI upper` = ci[2L],
               `Chance agreement` = pc, `Modified kappa` = kstar,
               Interpretation = if (!is.finite(kstar)) "" else if (kstar >= .74) "Excellent" else if (kstar >= .60) "Good" else if (kstar >= .40) "Fair" else "Poor",
               check.names = FALSE, stringsAsFactors = FALSE)
  })
  item <- do.call(rbind, rows)
  summary <- data.frame(
    Statistic = c("S-CVI/Ave", "S-CVI/UA", "Mean expert relevance score", "Rating cutoff", "Maximum rating"),
    Value = c(mean(item$`I-CVI`, na.rm = TRUE), mean(item$`I-CVI` == 1, na.rm = TRUE),
              mean(M, na.rm = TRUE), cutoff, max_score), stringsAsFactors = FALSE
  )
  list(item = item, summary = summary, data = M)
}

.r4vn_ts_known_groups <- function(scores, group, group_label, conf = 0.95) {
  g <- droplevels(as.factor(group))
  if (nlevels(g) < 2L) .r4vn_ts_stop("`known_groups` must contain at least two observed groups.")
  desc <- list(); tests <- list()
  for (sn in names(scores)) {
    x <- scores[[sn]]
    desc[[sn]] <- do.call(rbind, lapply(levels(g), function(lv) {
      z <- x[g == lv & is.finite(x)]
      data.frame(Scale = sn, Group = lv, N = length(z), Mean = mean(z), SD = stats::sd(z),
                 Median = stats::median(z), Q1 = stats::quantile(z, .25), Q3 = stats::quantile(z, .75),
                 check.names = FALSE, stringsAsFactors = FALSE)
    }))
    ok <- is.finite(x) & !is.na(g)
    if (nlevels(g) == 2L) {
      tt <- try(stats::t.test(x[ok] ~ g[ok], conf.level = conf), silent = TRUE)
      wt <- try(suppressWarnings(stats::wilcox.test(x[ok] ~ g[ok], exact = FALSE)), silent = TRUE)
      lev <- levels(g); a <- x[ok & g == lev[1L]]; b <- x[ok & g == lev[2L]]
      sp <- sqrt(((length(a)-1)*stats::var(a)+(length(b)-1)*stats::var(b))/(length(a)+length(b)-2))
      d <- (mean(b)-mean(a))/sp
      j <- 1 - 3/(4*(length(a)+length(b))-9)
      se_d <- sqrt((length(a)+length(b))/(length(a)*length(b)) + d^2/(2*(length(a)+length(b)-2)))
      qn <- stats::qnorm(1-(1-conf)/2)
      d_ci <- d+c(-1,1)*qn*se_d
      md_ci <- if (!inherits(tt,"try-error")) -rev(unname(tt$conf.int)) else c(NA_real_,NA_real_)
      tests[[sn]] <- data.frame(Scale = sn, Variable = group_label,
        `Group 1`=lev[1L], `Group 2`=lev[2L], Test = "Two-group comparison",
        `Parametric p` = if (!inherits(tt, "try-error")) tt$p.value else NA_real_,
        `Nonparametric p` = if (!inherits(wt, "try-error")) wt$p.value else NA_real_,
        `Mean difference (second-first)` = mean(b)-mean(a),
        `Mean difference CI lower`=md_ci[1L], `Mean difference CI upper`=md_ci[2L],
        `Cohen d` = d, `Cohen d CI lower`=d_ci[1L], `Cohen d CI upper`=d_ci[2L], `Hedges g` = j*d,
        `Eta squared` = NA_real_, check.names = FALSE, stringsAsFactors = FALSE)
    } else {
      fit <- stats::lm(x[ok] ~ g[ok]); av <- stats::anova(fit)
      kw <- stats::kruskal.test(x[ok] ~ g[ok])
      eta <- av$`Sum Sq`[1L]/sum(av$`Sum Sq`)
      tests[[sn]] <- data.frame(Scale = sn, Variable = group_label,
        `Group 1`=NA_character_, `Group 2`=NA_character_, Test = "Multiple-group comparison",
        `Parametric p` = av$`Pr(>F)`[1L], `Nonparametric p` = kw$p.value,
        `Mean difference (second-first)` = NA_real_,
        `Mean difference CI lower`=NA_real_, `Mean difference CI upper`=NA_real_,
        `Cohen d` = NA_real_, `Cohen d CI lower`=NA_real_, `Cohen d CI upper`=NA_real_, `Hedges g` = NA_real_,
        `Eta squared` = eta, check.names = FALSE, stringsAsFactors = FALSE)
    }
  }
  list(descriptive = do.call(rbind, desc), tests = do.call(rbind, tests), group = g, scores = scores)
}

.r4vn_ts_responsiveness <- function(before, after, conf = 0.95) {
  out <- list(); details <- list()
  for (sn in intersect(names(before), names(after))) {
    a <- before[[sn]]; b <- after[[sn]]
    z <- .r4vn_ts_agreement(a, b, sn, "Responsiveness", conf)
    d <- z$difference
    base_sd <- stats::sd(z$first)
    change_sd <- stats::sd(d)
    z$table$`Effect size` <- if (is.finite(base_sd) && base_sd > 0) mean(d)/base_sd else NA_real_
    z$table$SRM <- if (is.finite(change_sd) && change_sd > 0) mean(d)/change_sd else NA_real_
    out[[sn]] <- z$table
    details[[sn]] <- z
  }
  list(table = do.call(rbind, out), details = details)
}

.r4vn_ts_invariance <- function(X, original_data, factor_map, group_name,
                                ordered, estimator, missing, levels = c("configural", "metric", "scalar", "strict")) {
  if (is.null(group_name) || !nzchar(group_name)) .r4vn_ts_stop("Measurement invariance requires `cfa_group`.")
  if (!requireNamespace("lavaan", quietly = TRUE)) .r4vn_ts_stop("Measurement invariance requires the optional package `lavaan`.")
  levels <- match.arg(levels, several.ok = TRUE)
  levels <- c("configural", intersect(c("metric", "scalar", "strict"), levels))
  item_names <- colnames(X); model_items <- make.names(item_names, unique = TRUE)
  item_map <- stats::setNames(model_items, item_names)
  factor_names <- names(factor_map); model_factors <- make.names(factor_names, unique = TRUE)
  model <- paste(vapply(seq_along(factor_map), function(i) paste0(model_factors[i], " =~ ", paste(item_map[factor_map[[i]]], collapse = " + ")), character(1)), collapse = "\n")
  d <- as.data.frame(X, check.names = FALSE); names(d) <- model_items
  group_var <- ".r4vn_group"; d[[group_var]] <- original_data[[group_name]]
  ord <- if (isTRUE(ordered)) model_items else if (is.character(ordered)) unname(item_map[intersect(ordered, item_names)]) else character()
  est <- if (identical(estimator, "auto")) if (length(ord)) "WLSMV" else "MLR" else estimator
  scalar_equal <- if (length(ord)) c("loadings", "thresholds") else c("loadings", "intercepts")
  equality <- list(configural = NULL, metric = "loadings", scalar = scalar_equal,
                   strict = c(scalar_equal, "residuals"))
  fits <- list(); rows <- list()
  for (lv in levels) {
    args <- list(model = model, data = d, group = group_var, estimator = est, std.lv = TRUE)
    if (length(ord)) args$ordered <- ord else if (!is.null(missing) && nzchar(missing)) args$missing <- missing
    if (!is.null(equality[[lv]])) args$group.equal <- equality[[lv]]
    fit <- try(do.call(lavaan::cfa, args), silent = TRUE)
    if (inherits(fit, "try-error")) next
    fits[[lv]] <- fit
    fm <- lavaan::fitMeasures(fit, c("chisq", "df", "pvalue", "cfi", "tli", "rmsea", "srmr"))
    rows[[lv]] <- data.frame(Model = tools::toTitleCase(lv), ChiSquare = fm["chisq"], df = fm["df"], p = fm["pvalue"],
      CFI = fm["cfi"], TLI = fm["tli"], RMSEA = fm["rmsea"], SRMR = fm["srmr"], stringsAsFactors = FALSE)
  }
  if (!length(rows)) .r4vn_ts_stop("None of the requested measurement-invariance models converged.")
  tab <- do.call(rbind, rows); rownames(tab) <- NULL
  if (nrow(tab)) {
    tab$`Delta CFI` <- c(NA, diff(tab$CFI)); tab$`Delta RMSEA` <- c(NA, diff(tab$RMSEA)); tab$`Delta SRMR` <- c(NA, diff(tab$SRMR))
    tab$Supported <- ifelse(is.na(tab$`Delta CFI`), "Reference", ifelse(abs(tab$`Delta CFI`) <= .010 & tab$`Delta RMSEA` <= .015, "Yes", "No"))
  }
  list(table = tab, fits = fits, group = group_name, estimator = est)
}

.r4vn_ts_dif <- function(X, group, item_labels) {
  g <- droplevels(as.factor(group))
  if (nlevels(g) != 2L) .r4vn_ts_stop("Current DIF screening requires a two-level `cfa_group`.")
  total <- rowMeans(X, na.rm = TRUE)
  rows <- lapply(seq_along(X), function(j) {
    y <- X[[j]]; ok <- is.finite(y) & is.finite(total) & !is.na(g)
    if (sum(ok) < 20L) return(data.frame())
    fit0 <- stats::lm(y[ok] ~ total[ok])
    fit1 <- stats::lm(y[ok] ~ total[ok] + g[ok])
    fit2 <- stats::lm(y[ok] ~ total[ok] * g[ok])
    a01 <- stats::anova(fit0, fit1); a12 <- stats::anova(fit1, fit2)
    r20 <- summary(fit0)$r.squared; r21 <- summary(fit1)$r.squared; r22 <- summary(fit2)$r.squared
    data.frame(Item = item_labels[j], Variable = names(X)[j],
      `Uniform DIF p` = a01$`Pr(>F)`[2L], `Nonuniform DIF p` = a12$`Pr(>F)`[2L],
      `Uniform delta R2` = r21-r20, `Nonuniform delta R2` = r22-r21,
      Flag = ifelse(a01$`Pr(>F)`[2L] < .01 || a12$`Pr(>F)`[2L] < .01, "Review", "No signal"),
      check.names = FALSE, stringsAsFactors = FALSE)
  })
  do.call(rbind, rows)
}

.r4vn_ts_plot_draw <- function(spec) {
  old <- graphics::par(no.readonly = TRUE); on.exit(graphics::par(old), add = TRUE)
  # Use the platform sans-serif family consistently in the Plots pane and
  # raster output. This also preserves Vietnamese/Unicode labels where the
  # operating system provides the glyphs.
  graphics::par(family = "sans")
  tp <- spec$type
  if (tp == "item_diagnostics") {
    graphics::par(mar = c(5, max(7, min(16, max(nchar(spec$labels))*.55)), 3, 1))
    m <- rbind(spec$missing, spec$floor, spec$ceiling)
    graphics::barplot(m, beside = TRUE, horiz = TRUE, names.arg = spec$labels,
      col = c("#6BAED6", "#74C476", "#FD8D3C"), border = NA, las = 1,
      xlab = "Percent", main = spec$title)
    graphics::legend("bottomright", legend = c("Missing", "Floor", "Ceiling"),
      fill = c("#6BAED6", "#74C476", "#FD8D3C"), bty = "n")
  } else if (tp == "item_distribution") {
    graphics::par(mar = c(max(7, min(15, max(nchar(spec$labels))*.55)), 4, 3, 1))
    cols <- grDevices::hcl.colors(nrow(spec$proportions), "YlOrRd", rev = TRUE)
    graphics::barplot(spec$proportions, names.arg=spec$labels, las=2, col=cols,
      border=NA, ylab="Proportion", main=spec$title)
    graphics::legend("topright", legend=rownames(spec$proportions), fill=cols,
      title="Response", bty="n", cex=.8)
  } else if (tp == "reliability") {
    graphics::par(mar = c(max(7, min(15, max(nchar(spec$labels))*.55)), 4, 3, 1))
    ylim <- range(c(spec$item_total, spec$alpha_deleted, spec$alpha), finite = TRUE)
    graphics::plot(seq_along(spec$labels), spec$item_total, type = "b", pch = 16, col = "#2C7FB8",
      xaxt = "n", xlab = "", ylab = "Coefficient", ylim = ylim, main = spec$title)
    graphics::axis(1, at = seq_along(spec$labels), labels = spec$labels, las = 2, cex.axis = .75)
    graphics::lines(seq_along(spec$labels), spec$alpha_deleted, type = "b", pch = 17, col = "#D95F0E")
    graphics::abline(h = spec$alpha, lty = 2, col = "#555555")
    graphics::legend("bottomright", legend = c("Corrected item-total r", "Alpha if deleted", "Total alpha"),
      col = c("#2C7FB8", "#D95F0E", "#555555"), pch = c(16,17,NA), lty = c(1,1,2), bty = "n")
  } else if (tp == "correlation") {
    R <- spec$matrix; nr <- nrow(R)
    graphics::par(mar = c(max(7, min(15, max(nchar(spec$labels))*.55)), max(7, min(15, max(nchar(spec$labels))*.55)), 3, 2))
    graphics::image(seq_len(nr), seq_len(nr), t(R[nr:1, , drop=FALSE]), zlim=c(-1,1),
      col=grDevices::colorRampPalette(c("#B2182B","white","#2166AC"))(100), axes=FALSE, xlab="", ylab="", main=spec$title)
    graphics::axis(1, seq_len(nr), spec$labels, las=2, cex.axis=.7)
    graphics::axis(2, seq_len(nr), rev(spec$labels), las=2, cex.axis=.7)
  } else if (tp == "scores") {
    n <- length(spec$values); graphics::par(mfrow=c(ceiling(n/2), min(2,n)), mar=c(4,4,3,1))
    for (i in seq_along(spec$values)) graphics::hist(spec$values[[i]], col="#9ECAE1", border="white", main=names(spec$values)[i], xlab="Score")
  } else if (tp == "scree") {
    graphics::plot(spec$x, spec$observed, type="b", pch=16, col="#1F78B4", xlab="Component", ylab="Eigenvalue", main=spec$title)
    graphics::lines(spec$x, spec$random, type="b", pch=17, lty=2, col="#E31A1C")
    graphics::abline(h=1, lty=3, col="#666666"); graphics::legend("topright", c("Observed","Random 95th percentile"), col=c("#1F78B4","#E31A1C"), pch=c(16,17), lty=c(1,2), bty="n")
  } else if (tp == "loadings") {
    L <- spec$loadings; graphics::par(mar=c(max(7,min(15,max(nchar(rownames(L)))*.55)),4,3,1))
    graphics::matplot(seq_len(nrow(L)), L, type="b", pch=seq_len(ncol(L))+14, lty=1, xaxt="n", xlab="", ylab="Standardized loading", ylim=c(-1,1), main=spec$title)
    graphics::axis(1, seq_len(nrow(L)), rownames(L), las=2, cex.axis=.7); graphics::abline(h=0, col="#777777")
    graphics::legend("bottomright", colnames(L), col=seq_len(ncol(L)), pch=seq_len(ncol(L))+14, lty=1, bty="n")
  } else if (tp == "roc") {
    graphics::plot(c(0,1), c(0,1), type="n", xlab="1 - Specificity",
      ylab="Sensitivity", xlim=c(0,1), ylim=c(0,1), xaxs="i", yaxs="i",
      asp=1, main=spec$title)
    graphics::abline(0,1,lty=2,col="#888888")
    cols <- grDevices::hcl.colors(length(spec$curves), "Dark 3")
    for (i in seq_along(spec$curves)) {
      fpr <- pmax(0, pmin(1, 1-spec$curves[[i]]$specificity))
      tpr <- pmax(0, pmin(1, spec$curves[[i]]$sensitivity))
      graphics::lines(fpr, tpr, col=cols[i], lwd=2)
    }
    graphics::legend("bottomright", names(spec$curves), col=cols, lwd=2, bty="n")
  } else if (tp == "agreement") {
    graphics::par(mfrow=c(1,2), mar=c(4,4,3,1)); x<-spec$first; y<-spec$second
    graphics::plot(x,y,pch=16,col=grDevices::adjustcolor("#2C7FB8",.55),xlab="First measurement",ylab="Second measurement",main=paste(spec$title,"agreement")); graphics::abline(0,1,lty=2)
    avg<-(x+y)/2; d<-y-x; md<-mean(d); sdv<-stats::sd(d)
    graphics::plot(avg,d,pch=16,col=grDevices::adjustcolor("#D95F0E",.55),xlab="Mean of measurements",ylab="Difference",main="Bland-Altman plot"); graphics::abline(h=c(md,md-1.96*sdv,md+1.96*sdv),lty=c(1,2,2))
  } else if (tp == "known_groups") {
    graphics::boxplot(spec$score ~ spec$group, col="#A1D99B", xlab=spec$group_label, ylab="Score", main=spec$title)
  } else if (tp == "responsiveness") {
    graphics::plot(c(1,2), range(c(spec$first,spec$second),finite=TRUE), type="n", xaxt="n", xlab="", ylab="Score", main=spec$title)
    graphics::axis(1,c(1,2),c("Before","After")); graphics::segments(1,spec$first,2,spec$second,col=grDevices::adjustcolor("#3182BD",.2)); graphics::points(c(1,2),c(mean(spec$first),mean(spec$second)),pch=19,col="#CB181D",cex=1.4); graphics::lines(c(1,2),c(mean(spec$first),mean(spec$second)),col="#CB181D",lwd=3)
  } else if (tp == "coefficients") {
    graphics::par(mar=c(max(6,min(14,max(nchar(spec$labels))*.55)),4,3,1))
    graphics::barplot(spec$values, names.arg=spec$labels, las=2, ylim=c(min(0,min(spec$values,na.rm=TRUE)),1),
      col="#756BB1", border=NA, ylab="Coefficient", main=spec$title); graphics::abline(h=c(.5,.7,.9),lty=3,col="#999999")
  } else if (tp == "validity_forest") {
    y <- rev(seq_along(spec$estimate)); graphics::par(mar=c(4,max(8,min(18,max(nchar(spec$labels))*.5)),3,1))
    xr <- range(c(spec$lower,spec$upper,0),finite=TRUE); graphics::plot(spec$estimate,y,xlim=xr,ylim=c(.5,length(y)+.5),yaxt="n",pch=19,col="#238B45",xlab="Correlation (95% CI)",ylab="",main=spec$title)
    graphics::segments(spec$lower,y,spec$upper,y,col="#238B45",lwd=2); graphics::axis(2,y,spec$labels,las=2,cex.axis=.75); graphics::abline(v=0,lty=2,col="#777777")
  } else if (tp == "invariance") {
    d <- spec$data; keep <- is.finite(d$`Delta CFI`) | is.finite(d$`Delta RMSEA`) | is.finite(d$`Delta SRMR`)
    if (!any(keep)) { graphics::plot.new(); graphics::title(main=spec$title) } else {
      m <- rbind(abs(d$`Delta CFI`[keep]),abs(d$`Delta RMSEA`[keep]),abs(d$`Delta SRMR`[keep]))
      graphics::barplot(m,beside=TRUE,names.arg=d$Model[keep],col=c("#3182BD","#31A354","#E6550D"),border=NA,ylab="Absolute change",main=spec$title)
      graphics::abline(h=.01,lty=2,col="#555555"); graphics::legend("topleft",c("Delta CFI","Delta RMSEA","Delta SRMR"),fill=c("#3182BD","#31A354","#E6550D"),bty="n")
    }
  }
  invisible(spec)
}

.r4vn_ts_plot_svg <- function(spec, width = 8, height = 5.5) {
  path <- tempfile(fileext = ".svg")
  grDevices::svg(path, width = width, height = height, onefile = TRUE,
                 bg = "white", family = "sans")
  tryCatch(.r4vn_ts_plot_draw(spec), finally = grDevices::dev.off())
  txt <- paste(readLines(path, warn = FALSE, encoding = "UTF-8"), collapse = "\n")
  unlink(path)
  sub("^[\\s\\S]*?(<svg)", "\\1", txt, perl = TRUE)
}

.r4vn_ts_base64 <- function(x) {
  bytes <- as.integer(x)
  if (!length(bytes)) return("")
  alphabet <- strsplit(
    "ABCDEFGHIJKLMNOPQRSTUVWXYZabcdefghijklmnopqrstuvwxyz0123456789+/",
    "", fixed = TRUE
  )[[1L]]
  padding <- (3L - length(bytes) %% 3L) %% 3L
  if (padding) bytes <- c(bytes, rep.int(0L, padding))
  z <- matrix(bytes, ncol = 3L, byrow = TRUE)
  code <- cbind(
    bitwShiftR(z[, 1L], 2L),
    bitwOr(bitwShiftL(bitwAnd(z[, 1L], 3L), 4L), bitwShiftR(z[, 2L], 4L)),
    bitwOr(bitwShiftL(bitwAnd(z[, 2L], 15L), 2L), bitwShiftR(z[, 3L], 6L)),
    bitwAnd(z[, 3L], 63L)
  )
  encoded <- as.vector(t(matrix(alphabet[code + 1L], ncol = 4L)))
  if (padding) encoded[(length(encoded) - padding + 1L):length(encoded)] <- "="
  paste0(encoded, collapse = "")
}

.r4vn_ts_plot_png <- function(spec, width = 1152L, height = 792L,
                               res = 144L) {
  path <- tempfile(fileext = ".png")
  on.exit(unlink(path), add = TRUE)
  args <- list(filename = path, width = width, height = height, units = "px",
               res = res, bg = "white")
  if (isTRUE(capabilities("cairo"))) args$type <- "cairo-png"
  do.call(grDevices::png, args)
  tryCatch(.r4vn_ts_plot_draw(spec), finally = grDevices::dev.off())
  size <- file.info(path)$size
  if (!is.finite(size) || size <= 0) .r4vn_ts_stop("Could not render the Viewer plot as PNG.")
  bytes <- readBin(path, what = "raw", n = size)
  paste0(
    "<img class=\"ts-plot-image\" alt=\"",
    .r4vn_ts_escape(.r4vn_ts_or(spec$title, "tabscale plot")),
    "\" src=\"data:image/png;base64,", .r4vn_ts_base64(bytes), "\">"
  )
}

#' Comprehensive Scale Analysis
#'
#' @description
#' `tabscale()` provides a one-command, publication-ready psychometric report.
#' It covers internal consistency, stability, equivalence, inter-rater
#' reliability, measurement error, content validity, structural validity,
#' convergent/discriminant and known-groups validity, criterion validity,
#' measurement invariance, DIF screening, and responsiveness when the required
#' data are supplied. Variable labels are used throughout whenever available.
#'
#' @usage
#' tabscale(data = NULL, vars = NULL, factor = NULL, reverse = NULL,
#'          report = c("auto", "brief", "full", "custom"),
#'          range = NULL, score = c("mean", "sum"), min_valid = NULL,
#'          missing = c("pairwise", "complete"), cor_method = c("pearson", "spearman"),
#'          reliability = TRUE, bootstrap = 0, conf = 0.95,
#'          retest = NULL, parallel_form = NULL, raters = NULL,
#'          rater_type = c("auto", "continuous", "categorical"),
#'          content = NULL, content_cutoff = 3, content_max = 4,
#'          validity = TRUE,
#'          efa = NULL, efa_method = c("pa", "ml"), nfactor = NULL,
#'          rotation = c("varimax", "promax", "none"), parallel_iter = 100,
#'          cfa = NULL, ordered = FALSE, estimator = "auto",
#'          cfa_missing = "fiml", cfa_group = NULL,
#'          modification = FALSE, modification_min = 10,
#'          invariance = FALSE,
#'          invariance_levels = c("configural", "metric", "scalar", "strict"),
#'          dif = FALSE, convergent = NULL, discriminant = NULL,
#'          convergent_min = 0.50, discriminant_max = 0.30,
#'          known_groups = NULL,
#'          gold = NULL, event = NULL, direction = c("auto", "higher", "lower"),
#'          post = NULL,
#'          name = FALSE, digit = 2, p_digit = 3,
#'          template = c("journal", "clean", "minimal"),
#'          plot = TRUE, plot_types = "auto",
#'          viewer_plot_format = c("png", "svg"),
#'          append = NULL, file = NULL,
#'          title = NULL, interpretation = FALSE,
#'          raw = TRUE, show = TRUE, seed = NULL)
#'
#' @param data Optional data frame. When omitted, the active R4VN data frame is used.
#' @param vars Items created by `vars()`, a character vector, or a one-sided formula.
#' @param factor Optional named list defining subscales/CFA factors.
#' @param reverse Optional items to reverse-score.
#' @param report Output profile. `"auto"` runs reliability plus EFA and any
#'   data-dependent modules requested by their arguments; `"brief"` omits
#'   factor models; `"full"` also runs CFA when a factor map and lavaan are
#'   available; `"custom"` follows the module switches exactly.
#' @param range Two numeric values giving the minimum and maximum item score.
#' @param score Calculate scale scores as the item `"mean"` or `"sum"`.
#' @param min_valid Minimum valid items. A value in `(0,1]` is treated as a proportion.
#'   A named vector can define separate minima for factors and `Total`.
#' @param missing Correlation/covariance handling: pairwise or complete observations.
#' @param cor_method Pearson or Spearman item correlations.
#' @param reliability Logical; calculate reliability statistics.
#' @param bootstrap Number of nonparametric bootstrap replicates for alpha and
#'   omega-total confidence intervals. Zero uses a Feldt interval for alpha.
#' @param conf Confidence level for alpha and AUC intervals.
#' @param retest Items measured again, in the same order as `vars()`, for
#'   test-retest reliability, ICC, SEM, MDC, and Bland-Altman analysis.
#' @param parallel_form Items from an equivalent form, in the same order as
#'   `vars()`, for parallel-form reliability.
#' @param raters Two or more variables containing ratings of the same subjects.
#' @param rater_type Treat ratings as continuous or categorical; `"auto"`
#'   chooses categorical for variables with at most 10 observed levels.
#' @param content Expert-by-item matrix/data frame of content-relevance ratings.
#' @param content_cutoff Minimum rating counted as content-relevant.
#' @param content_max Maximum possible content rating, retained in the report.
#' @param validity Logical master switch for validity modules.
#' @param efa Logical or `NULL`; run EFA, KMO, Bartlett, and parallel analysis.
#' @param efa_method Principal-axis (`"pa"`) or maximum-likelihood (`"ml"`) EFA.
#' @param nfactor Number of EFA factors. When `NULL`, parallel analysis is used.
#' @param rotation EFA rotation.
#' @param parallel_iter Number of Monte Carlo samples for parallel analysis.
#' @param cfa Logical or `NULL`; run CFA using the optional `lavaan` package.
#' @param ordered Logical or character item names treated as ordinal in CFA.
#' @param estimator CFA estimator. `"auto"` uses WLSMV for ordered items and MLR otherwise.
#' @param cfa_missing Missing-data option passed to lavaan for non-ordinal CFA.
#' @param cfa_group Optional grouping variable name for multiple-group CFA.
#' @param modification Logical; include large CFA modification indices.
#' @param modification_min Minimum modification index displayed.
#' @param invariance Logical; test configural, metric, scalar, and strict
#'   measurement invariance across `cfa_group`.
#' @param invariance_levels Invariance levels to fit.
#' @param dif Logical; screen uniform and non-uniform differential item
#'   functioning across a two-level `cfa_group`.
#' @param convergent External variables used for convergent validity.
#' @param discriminant External variables used for discriminant validity.
#' @param convergent_min Prespecified minimum absolute convergent correlation.
#' @param discriminant_max Prespecified maximum absolute discriminant correlation.
#' @param known_groups Grouping variable for known-groups validity.
#' @param gold Optional criterion or binary gold-standard variable.
#' @param event Event level for binary gold-standard ROC analysis.
#' @param direction Whether higher or lower scores predict the event; `"auto"` chooses the direction with AUC at least 0.5.
#' @param post Post-intervention/follow-up items, in the same order as `vars()`,
#'   for responsiveness (effect size and standardized response mean).
#' @param name `FALSE` does not modify data. `TRUE` creates `scale_total` and
#'   subscale variables; a character value supplies the score prefix.
#' @param digit Decimal places for estimates.
#' @param p_digit Decimal places for p-values.
#' @param template HTML style.
#' @param plot Logical; create all applicable graphics in both the R Plots pane
#'   and the HTML Viewer.
#' @param plot_types `"auto"`, `"all"`, or any of `"items"`, `"distributions"`,
#'   `"correlation"`, `"reliability"`, `"scores"`, `"scree"`,
#'   `"loadings"`, `"cfa_loadings"`, `"roc"`, `"test_retest"`,
#'   `"parallel_form"`, `"inter_rater"`, `"content"`,
#'   `"external_validity"`, `"known_groups"`, `"invariance"`, and
#'   `"responsiveness"`. The automatic set omits secondary plots that often
#'   add visual clutter: item completeness, item-response distributions,
#'   item-reliability diagnostics, external-validity correlations,
#'   known-groups boxplots, and individual responsiveness trajectories.
#'   Their statistical tables remain in the report. Request any of these names
#'   explicitly, or use `plot_types = "all"`, when the graphic is needed.
#' @param viewer_plot_format Format used to embed plots in the HTML Viewer.
#'   The default `"png"` uses high-resolution, self-contained images and gives
#'   consistent fonts in RStudio Viewer, Chrome, and saved HTML files. `"svg"`
#'   retains vector graphics but may render text differently across browsers
#'   because SVG font substitution is controlled by the local system.
#' @param append Optional previous R4VN table object or HTML file.
#' @param file Optional HTML output path.
#' @param title Optional table title.
#' @param interpretation Logical; add cautious automatic interpretation. The
#'   default is `FALSE`.
#' @param raw Logical; retain numerical result components.
#' @param show Logical; open the HTML result.
#' @param seed Optional random seed used by parallel analysis. The default `NULL` does not set a seed.
#'
#' @details
#' The default `report = "auto"` produces descriptive item distributions,
#' missing/floor/ceiling effects, corrected item-total correlations, alpha with
#' 95% CI, standardized alpha, omega total, split-half coefficients, all six
#' Guttman lambdas, KMO, Bartlett's test, parallel analysis, EFA, score
#' distributions, and matching graphics. KR-20 is added for binary items.
#' Ordinal alpha, omega hierarchical, and the greatest lower bound are added
#' when the optional package `psych` is installed.
#'
#' Additional data activate stability/test-retest, parallel-form, inter-rater,
#' content, convergent, discriminant, known-groups, criterion, responsiveness,
#' measurement-invariance, and DIF sections. CFA and invariance use `lavaan`.
#' Face validity is inherently qualitative and is therefore identified in the
#' coverage table rather than assigned a spurious numeric coefficient.
#'
#' The same plot specifications are rendered in the R graphics device and in
#' the HTML Viewer. In RStudio, use the Plots pane arrows to review every graph,
#' or rerun selected graphs with `plot(result, which = "roc")`.
#'
#' @return Invisibly returns an object inheriting from `r4vn_tabscale` and
#' `r4vn_tab`. Components include `tables`, `coverage`, `descriptive`,
#' `reliability`, `scores`, `test_retest`, `parallel_form`, `inter_rater`,
#' `content_validity`, `efa`, `cfa`, `convergent_validity`,
#' `discriminant_validity`, `known_groups`, `criterion_validity`, `invariance`,
#' `dif`, `responsiveness`, `plots`, `interpretation`, and `file`.
#'
#' @seealso `vars`, `tab`, `tabmulti`, `tabexport`
#' @family R4VN tables
#'
#' @examples
#' set.seed(2026)
#' n <- 120
#' f1 <- rnorm(n)
#' f2 <- 0.35 * f1 + rnorm(n, sd = 0.94)
#' make_item <- function(z) as.integer(cut(z, quantile(z, 0:5/5),
#'                                        include.lowest = TRUE, labels = FALSE))
#' latent <- list(0.8*f1, 0.7*f1, 0.9*f1, -0.7*f1,
#'                0.8*f2, 0.7*f2, 0.9*f2, 0.6*f2)
#' base_items <- lapply(latent, function(z) make_item(z + rnorm(n)))
#' dat <- as.data.frame(base_items)
#' names(dat) <- paste0("q", 1:8)
#' for (j in 1:8) {
#'   dat[[paste0("q", j, "_retest")]] <- pmax(
#'     1,
#'     pmin(
#'       5,
#'       dat[[paste0("q", j)]] + sample(-1:1, n, TRUE, c(.1, .8, .1))
#'     )
#'   )
#'   dat[[paste0("q", j, "_formb")]] <- make_item(latent[[j]] + rnorm(n))
#'   dat[[paste0("q", j, "_post")]] <- pmax(1, pmin(5, dat[[paste0("q",j)]] + rbinom(n,1,.35)))
#'   attr(dat[[paste0("q",j)]], "label") <- paste("Well-being item", j)
#' }
#' dat$convergent_measure <- f1 + f2 + rnorm(n, sd=.6)
#' dat$unrelated_measure <- rnorm(n)
#' dat$known_group <- factor(ifelse(f1+f2>0,"Higher expected score","Lower expected score"))
#' dat$gold <- factor(ifelse(f1 + f2 + rnorm(n) > 0, "Yes", "No"),
#'                    levels = c("No", "Yes"))
#' dat$rater1 <- sample(1:4,n,TRUE); dat$rater2 <- dat$rater1
#' dat$rater3 <- dat$rater1
#' dat$rater2[sample(n,30)] <- sample(1:4,30,TRUE)
#' dat$rater3[sample(n,35)] <- sample(1:4,35,TRUE)
#' attr(dat$known_group,"label") <- "Prespecified clinical group"
#' attr(dat$gold,"label") <- "Clinical gold standard"
#'
#' # 1. One-command automatic report: reliability, factorability, EFA, plots.
#' tb <- tabscale(
#'   dat,
#'   vars = vars(q1, q2, q3, q4, q5, q6, q7, q8),
#'   factor = list(Domain1 = vars(q1, q2, q3, q4),
#'                 Domain2 = vars(q5, q6, q7, q8)),
#'   reverse = vars(q4), range = c(1, 5),
#'   nfactor = 2, parallel_iter = 10,
#'   plot = FALSE, show = FALSE
#' )
#' tb$reliability_summary
#' tb$efa$loadings
#' if (interactive()) {
#'   plot(tb)                         # all plots in the Plots pane
#'   plot(tb, which = "reliability") # one selected plot
#' }
#'
#' # 2. Stability, parallel forms, inter-rater reliability, and measurement error.
#' rel <- tabscale(
#'   dat, vars=vars(q1,q2,q3,q4,q5,q6,q7,q8), reverse=vars(q4), range=c(1,5),
#'   retest=vars(q1_retest,q2_retest,q3_retest,q4_retest,
#'               q5_retest,q6_retest,q7_retest,q8_retest),
#'   parallel_form=vars(q1_formb,q2_formb,q3_formb,q4_formb,
#'                      q5_formb,q6_formb,q7_formb,q8_formb),
#'   raters=vars(rater1,rater2,rater3),
#'   plot=FALSE, show=FALSE)
#' rel$test_retest$table
#' rel$parallel_form$table
#' rel$inter_rater$table
#'
#' # 3. Convergent, discriminant, known-groups, criterion validity,
#' #    responsiveness, and automatic ROC curves.
#' val <- tabscale(
#'   dat, vars=vars(q1,q2,q3,q4,q5,q6,q7,q8), reverse=vars(q4), range=c(1,5),
#'   convergent=vars(convergent_measure), discriminant=vars(unrelated_measure),
#'   known_groups=known_group, gold=gold, event="Yes",
#'   post=vars(q1_post,q2_post,q3_post,q4_post,q5_post,q6_post,q7_post,q8_post),
#'   interpretation=TRUE, plot=FALSE, show=FALSE)
#' val$convergent_validity
#' val$discriminant_validity
#' val$known_groups$tests
#' val$criterion_validity$table
#' val$responsiveness$table
#'
#' # 4. Content validity: experts in rows and items in columns.
#' expert_ratings <- as.data.frame(matrix(sample(2:4, 6*8, TRUE,
#'   prob=c(.10,.30,.60)), nrow=6, dimnames=list(NULL,paste0("q",1:8))))
#' content_result <- tabscale(
#'   dat, vars=vars(q1,q2,q3,q4,q5,q6,q7,q8), range=c(1,5),
#'   content=expert_ratings, content_cutoff=3,
#'   plot=FALSE, show=FALSE)
#' content_result$content_validity$summary
#' content_result$content_validity$item
#'
#' # 5. CFA, composite reliability, AVE, HTMT, Fornell-Larcker,
#' #    measurement invariance, and DIF screening.
#' if (interactive() && requireNamespace("lavaan", quietly = TRUE)) {
#'   cfa_result <- tabscale(
#'     dat, vars=vars(q1,q2,q3,q4,q5,q6,q7,q8),
#'     factor=list(Domain1=vars(q1,q2,q3,q4), Domain2=vars(q5,q6,q7,q8)),
#'     reverse=vars(q4), range=c(1,5), cfa=TRUE, ordered=TRUE,
#'     cfa_group=known_group, invariance=TRUE, dif=TRUE,
#'     modification=TRUE, plot=FALSE, show=FALSE)
#'   cfa_result$cfa$reliability
#'   cfa_result$cfa$htmt
#'   cfa_result$invariance$table
#'   cfa_result$dif
#' }
#' @export
tabscale <- function(data = NULL, vars = NULL, factor = NULL, reverse = NULL,
                      report = c("auto", "brief", "full", "custom"),
                      range = NULL, score = c("mean", "sum"), min_valid = NULL,
                      missing = c("pairwise", "complete"), cor_method = c("pearson", "spearman"),
                      reliability = TRUE, bootstrap = 0L, conf = 0.95,
                      retest = NULL, parallel_form = NULL, raters = NULL,
                      rater_type = c("auto", "continuous", "categorical"),
                      content = NULL, content_cutoff = 3, content_max = 4,
                      validity = TRUE,
                      efa = NULL, efa_method = c("pa", "ml"), nfactor = NULL,
                      rotation = c("varimax", "promax", "none"), parallel_iter = 100L,
                      cfa = NULL, ordered = FALSE, estimator = "auto",
                      cfa_missing = "fiml", cfa_group = NULL,
                      modification = FALSE, modification_min = 10,
                      invariance = FALSE,
                      invariance_levels = c("configural", "metric", "scalar", "strict"),
                      dif = FALSE,
                      convergent = NULL, discriminant = NULL,
                      convergent_min = 0.50, discriminant_max = 0.30,
                      known_groups = NULL,
                      gold = NULL, event = NULL, direction = c("auto", "higher", "lower"),
                      post = NULL,
                      name = FALSE, digit = 2L, p_digit = 3L,
                      template = c("journal", "clean", "minimal"),
                      plot = TRUE, plot_types = "auto",
                      viewer_plot_format = c("png", "svg"),
                      append = NULL, file = NULL, title = NULL,
                      interpretation = FALSE, raw = TRUE, show = TRUE, seed = NULL) {
  env <- parent.frame()
  data_supplied <- !missing(data)
  data_expr <- substitute(data)
  gold_expr <- substitute(gold)
  cfa_group_expr <- substitute(cfa_group)
  known_groups_expr <- substitute(known_groups)

  # Allow tabscale(vars(...)) in active-data mode while preserving
  # tabscale(data, vars = vars(...)) for explicit data.
  if (!is.null(data) && !is.data.frame(data) && is.null(vars) && (inherits(data, "r4vn_vars") || is.character(data) || inherits(data, "formula"))) {
    vars <- data
    data <- NULL
    data_expr <- quote(NULL)
    data_supplied <- FALSE
  }

  data <- .r4vn_ts_get_data(data)
  report <- match.arg(report)
  score <- match.arg(score)
  missing <- match.arg(missing)
  cor_method <- match.arg(cor_method)
  efa_method <- match.arg(efa_method)
  rotation <- match.arg(rotation)
  direction <- match.arg(direction)
  rater_type <- match.arg(rater_type)
  template <- match.arg(template)
  viewer_plot_format <- match.arg(viewer_plot_format)
  valid_plot_types <- c("items", "distributions", "correlation", "reliability", "scores", "scree",
    "loadings", "cfa_loadings", "roc", "test_retest", "parallel_form", "inter_rater",
    "content", "external_validity", "known_groups", "invariance", "responsiveness")
  if (!identical(plot_types, "auto") && !identical(plot_types, "all")) {
    plot_types <- unique(as.character(plot_types))
    bad_plot_types <- setdiff(plot_types, valid_plot_types)
    if (length(bad_plot_types)) .r4vn_ts_stop("Unknown `plot_types`: ", paste(bad_plot_types, collapse = ", "))
  }
  digit <- max(0L, as.integer(digit))
  p_digit <- max(0L, as.integer(p_digit))
  if (!is.finite(conf) || conf <= 0 || conf >= 1) .r4vn_ts_stop("`conf` must be between 0 and 1.")
  if (is.null(efa)) efa <- report %in% c("auto", "full")
  if (is.null(cfa)) {
    wants_auto_cfa <- isTRUE(report == "full") && !is.null(factor)
    cfa <- wants_auto_cfa && requireNamespace("lavaan", quietly = TRUE)
    if (wants_auto_cfa && !isTRUE(cfa)) .r4vn_ts_warn("`report = \"full\"` could not run CFA because the optional package `lavaan` is not installed.")
  }
  if (report == "brief") efa <- cfa <- FALSE
  validity_enabled <- isTRUE(validity)

  factor_map0 <- .r4vn_ts_factor_map(factor)
  item_names <- if (is.null(vars)) unique(unlist(factor_map0, use.names = FALSE)) else .r4vn_ts_names(vars, "vars")
  if (!length(item_names)) .r4vn_ts_stop("No scale items were supplied.")
  if (anyDuplicated(item_names)) item_names <- unique(item_names)
  bad <- setdiff(item_names, names(data))
  if (length(bad)) .r4vn_ts_stop("Items not found in data: ", paste(bad, collapse = ", "))
  if (length(item_names) < 2L) .r4vn_ts_stop("At least two scale items are required.")
  factor_map <- .r4vn_ts_factor_map(factor, item_names)

  reverse_names <- .r4vn_ts_names(reverse, "reverse", allow_null = TRUE)
  bad_reverse <- setdiff(reverse_names, item_names)
  if (length(bad_reverse)) .r4vn_ts_stop("Reverse-coded items not found in `vars`: ", paste(bad_reverse, collapse = ", "))
  if (length(reverse_names) && (is.null(range) || length(range) != 2L || any(!is.finite(range)))) {
    .r4vn_ts_stop("Supply `range = c(minimum, maximum)` when reverse-scored items are used.")
  }
  if (!is.null(range)) {
    range <- sort(as.numeric(range))
    if (length(range) != 2L || any(!is.finite(range)) || range[1L] >= range[2L]) .r4vn_ts_stop("`range` must contain two increasing finite values.")
  }

  converted <- lapply(item_names, function(v) .r4vn_ts_numeric_item(data[[v]], v))
  names(converted) <- item_names
  X_raw <- as.data.frame(stats::setNames(lapply(converted, `[[`, "value"), item_names), check.names = FALSE)
  X <- X_raw

  if (!is.null(range)) {
    outside <- vapply(X, function(z) any(!is.na(z) & (z < range[1L] | z > range[2L])), logical(1))
    if (any(outside)) .r4vn_ts_warn("Values outside `range` were set to missing for: ", paste(names(outside)[outside], collapse = ", "))
    for (v in names(outside)[outside]) {
      X_raw[[v]][!is.na(X_raw[[v]]) & (X_raw[[v]] < range[1L] | X_raw[[v]] > range[2L])] <- NA_real_
      X[[v]][!is.na(X[[v]]) & (X[[v]] < range[1L] | X[[v]] > range[2L])] <- NA_real_
    }
  }

  if (length(reverse_names)) {
    lo <- range[1L]
    hi <- range[2L]
    for (v in reverse_names) X[[v]] <- lo + hi - X[[v]]
  }

  variances <- vapply(X, stats::var, numeric(1), na.rm = TRUE)
  no_variation <- names(variances)[!is.finite(variances) | variances <= 0]
  if (length(no_variation)) .r4vn_ts_stop("Items with no usable variation: ", paste(no_variation, collapse = ", "))

  # Item descriptive table --------------------------------------------------
  response_values <- sort(unique(unlist(X_raw, use.names = FALSE)))
  response_values <- response_values[is.finite(response_values)]
  if (length(response_values) > 15L) .r4vn_ts_warn("More than 15 distinct response values were found; the item distribution table will be wide.")
  response_columns <- vapply(response_values, function(value) {
    labs <- unique(unlist(lapply(converted, function(z) {
      if (is.null(z$labels) || !as.character(value) %in% names(z$labels)) return(character())
      as.character(z$labels[[as.character(value)]])
    }), use.names = FALSE))
    labs <- labs[!is.na(labs) & nzchar(labs)]
    if (length(labs) == 1L && labs != as.character(value)) paste0(value, ": ", labs) else paste0("Score ", value)
  }, character(1))
  response_columns <- make.unique(response_columns)
  descriptive_rows <- lapply(seq_along(item_names), function(i) {
    v <- item_names[i]
    x0 <- X_raw[[v]]
    xs <- X[[v]]
    n <- sum(!is.na(xs))
    miss <- sum(is.na(xs))
    q <- if (n) stats::quantile(xs, c(.25, .5, .75), na.rm = TRUE, names = FALSE) else rep(NA_real_, 3L)
    minv <- if (n) min(xs, na.rm = TRUE) else NA_real_
    maxv <- if (n) max(xs, na.rm = TRUE) else NA_real_
    floor_value <- if (!is.null(range)) range[1L] else minv
    ceiling_value <- if (!is.null(range)) range[2L] else maxv
    base <- data.frame(
      Item = .r4vn_ts_label(data[[v]], v),
      Variable = v,
      Subscale = {
        z <- names(factor_map)[vapply(factor_map, function(set) v %in% set, logical(1))]
        if (length(z)) paste(z, collapse = ", ") else ""
      },
      Reversed = if (v %in% reverse_names) "Yes" else "No",
      N = n,
      Missing = sprintf("%d (%.1f%%)", miss, 100 * miss / nrow(data)),
      `Mean (SD)` = if (n) sprintf(paste0("%.", digit, "f (%.", digit, "f)"), mean(xs, na.rm = TRUE), stats::sd(xs, na.rm = TRUE)) else "",
      `Median (IQR)` = if (n) sprintf(paste0("%.", digit, "f (%.", digit, "f-%.", digit, "f)"), q[2L], q[1L], q[3L]) else "",
      `Range` = if (n) sprintf(paste0("%.", digit, "f-%.", digit, "f"), minv, maxv) else "",
      `Floor %` = if (n) 100 * mean(xs == floor_value, na.rm = TRUE) else NA_real_,
      `Ceiling %` = if (n) 100 * mean(xs == ceiling_value, na.rm = TRUE) else NA_real_,
      check.names = FALSE,
      stringsAsFactors = FALSE
    )
    for (j in seq_along(response_values)) {
      val <- response_values[j]
      nn <- sum(x0 == val, na.rm = TRUE)
      den <- sum(!is.na(x0))
      base[[response_columns[j]]] <- if (den) sprintf("%d (%.1f%%)", nn, 100 * nn / den) else ""
    }
    base
  })
  descriptive_table <- do.call(rbind, descriptive_rows)

  # Scores and reliability --------------------------------------------------
  score_objects <- .r4vn_ts_subscale_scores(X, factor_map, score, min_valid)
  score_df <- as.data.frame(lapply(score_objects, `[[`, "score"), check.names = FALSE)
  names(score_df) <- names(score_objects)
  item_labels <- vapply(item_names, function(v) .r4vn_ts_label(data[[v]], v), character(1))

  alternate_scores <- function(specification, argument) {
    nm <- .r4vn_ts_names(specification, argument, allow_null = TRUE)
    if (!length(nm)) return(NULL)
    bad <- setdiff(nm, names(data))
    if (length(bad)) .r4vn_ts_stop("Variables in `", argument, "` not found in data: ", paste(bad, collapse = ", "))
    if (length(nm) != length(item_names)) {
      .r4vn_ts_stop("`", argument, "` must contain ", length(item_names), " items in the same order as `vars`.")
    }
    conv <- lapply(nm, function(v) .r4vn_ts_numeric_item(data[[v]], v)$value)
    Z <- as.data.frame(stats::setNames(conv, item_names), check.names = FALSE)
    if (!is.null(range)) {
      for (v in item_names) Z[[v]][!is.na(Z[[v]]) & (Z[[v]] < range[1L] | Z[[v]] > range[2L])] <- NA_real_
    }
    if (length(reverse_names)) {
      for (v in reverse_names) Z[[v]] <- range[1L] + range[2L] - Z[[v]]
    }
    obj <- .r4vn_ts_subscale_scores(Z, factor_map, score, min_valid)
    list(items = nm, matrix = Z, scores = as.data.frame(lapply(obj, `[[`, "score"), check.names = FALSE))
  }

  retest_result <- alternate_scores(retest, "retest")
  parallel_result <- alternate_scores(parallel_form, "parallel_form")
  post_result <- if (validity_enabled) alternate_scores(post, "post") else NULL

  stability <- equivalence <- responsiveness <- NULL
  if (!is.null(retest_result)) {
    details <- lapply(names(score_df), function(nm) .r4vn_ts_agreement(score_df[[nm]], retest_result$scores[[nm]], nm, "Test-retest reliability", conf))
    names(details) <- names(score_df)
    stability <- list(table = do.call(rbind, lapply(details, `[[`, "table")), details = details,
                      items = retest_result$items, scores = retest_result$scores)
  }
  if (!is.null(parallel_result)) {
    details <- lapply(names(score_df), function(nm) .r4vn_ts_agreement(score_df[[nm]], parallel_result$scores[[nm]], nm, "Parallel-form reliability", conf))
    names(details) <- names(score_df)
    equivalence <- list(table = do.call(rbind, lapply(details, `[[`, "table")), details = details,
                        items = parallel_result$items, scores = parallel_result$scores)
  }
  if (!is.null(post_result)) responsiveness <- .r4vn_ts_responsiveness(score_df, post_result$scores, conf)

  content_result <- if (validity_enabled) .r4vn_ts_content_validity(content, content_cutoff, content_max, conf) else NULL

  rater_result <- NULL
  rater_names <- .r4vn_ts_names(raters, "raters", allow_null = TRUE)
  if (length(rater_names)) {
    bad <- setdiff(rater_names, names(data))
    if (length(bad)) .r4vn_ts_stop("Rater variables not found in data: ", paste(bad, collapse = ", "))
    if (length(rater_names) < 2L) .r4vn_ts_stop("`raters` must contain at least two rating variables.")
    rater_data <- data[, rater_names, drop = FALSE]
    numeric_raters <- as.data.frame(lapply(rater_data, function(z) .r4vn_ts_numeric_item(z, "rater")$value), check.names = FALSE)
    names(numeric_raters) <- rater_names
    type_use <- rater_type
    if (type_use == "auto") type_use <- if (all(vapply(rater_data, function(z) length(unique(z[!is.na(z)])) <= 10L, logical(1)))) "categorical" else "continuous"
    icc <- .r4vn_ts_icc(numeric_raters, conf)
    kap <- if (type_use == "categorical") .r4vn_ts_kappa(rater_data, weighted = FALSE, conf = conf) else
      c(N=unname(icc["N"]), Raters=length(rater_names), Agreement=NA, Kappa=NA, lower=NA, upper=NA)
    kap_weighted <- if (type_use == "categorical" && length(rater_names) == 2L)
      .r4vn_ts_kappa(rater_data, weighted = TRUE, conf = conf) else kap
    rater_table <- data.frame(
      Type = type_use, N = unname(icc["N"]), Raters = length(rater_names),
      `Exact agreement` = unname(kap["Agreement"]),
      `Unweighted/Fleiss kappa` = unname(kap["Kappa"]),
      `Kappa CI lower` = unname(kap["lower"]), `Kappa CI upper` = unname(kap["upper"]),
      `Weighted kappa (two ordinal raters)` = unname(kap_weighted["Kappa"]),
      `Weighted kappa CI lower` = unname(kap_weighted["lower"]),
      `Weighted kappa CI upper` = unname(kap_weighted["upper"]),
      `ICC absolute single` = unname(icc["ICC_A1"]), `ICC absolute average` = unname(icc["ICC_Ak"]),
      `ICC consistency single` = unname(icc["ICC_C1"]), `ICC consistency average` = unname(icc["ICC_Ck"]),
      `ICC CI lower` = unname(icc["lower"]), `ICC CI upper` = unname(icc["upper"]),
      SEM = unname(icc["SEM"]), MDC95 = unname(icc["MDC95"]),
      check.names = FALSE, stringsAsFactors = FALSE
    )
    rater_result <- list(type = type_use, variables = rater_names,
                         labels = vapply(rater_names, function(v) .r4vn_ts_label(data[[v]], v), character(1)),
                         table = rater_table, data = rater_data, icc = icc,
                         kappa = kap, weighted_kappa = kap_weighted)
  }

  reliability_results <- list()
  reliability_summary <- data.frame()
  item_analysis <- data.frame()
  if (isTRUE(reliability)) {
    scale_sets <- c(factor_map, list(Total = item_names))
    reliability_results <- lapply(scale_sets, function(set) .r4vn_ts_alpha(X[, set, drop = FALSE], missing, cor_method, bootstrap, conf))
    reliability_summary <- do.call(rbind, lapply(names(scale_sets), function(nm) {
      set <- scale_sets[[nm]]
      rr <- reliability_results[[nm]]
      sc <- score_objects[[nm]]$score
      data.frame(
        Scale = nm,
        Items = length(set),
        ValidScores = sum(is.finite(sc)),
        `Score mean` = .r4vn_ts_safe(sc, mean),
        `Score SD` = .r4vn_ts_safe(sc, stats::sd),
        `Score min` = .r4vn_ts_safe(sc, min),
        `Score max` = .r4vn_ts_safe(sc, max),
        Alpha = rr$alpha,
        `Alpha CI lower` = rr$alpha_ci[1L],
        `Alpha CI upper` = rr$alpha_ci[2L],
        `Alpha CI method` = rr$alpha_ci_method,
        `Standardized alpha` = rr$alpha_std,
        `Ordinal alpha` = rr$ordinal_alpha,
        `Omega total` = rr$omega_total,
        `Omega CI lower` = rr$omega_ci[1L],
        `Omega CI upper` = rr$omega_ci[2L],
        `Omega hierarchical` = rr$omega_hierarchical,
        `Greatest lower bound` = rr$glb,
        KR20 = rr$kr20,
        `Mean inter-item r` = rr$mean_interitem,
        `Median inter-item r` = rr$median_interitem,
        `Split-half r` = rr$split_half_r,
        `Spearman-Brown` = rr$spearman_brown,
        `Mean split-half reliability` = rr$split_half_mean,
        `Best split-half reliability (lambda 4)` = rr$split_half_best,
        `Guttman lambda 1` = rr$lambda1,
        `Guttman lambda 2` = rr$lambda2,
        `Guttman lambda 3` = rr$lambda3,
        `Guttman lambda 4` = rr$lambda4,
        `Guttman lambda 5` = rr$lambda5,
        `Guttman lambda 6` = rr$lambda6,
        `Minimum valid items` = score_objects[[nm]]$min_valid,
        check.names = FALSE,
        stringsAsFactors = FALSE
      )
    }))

    rr_total <- reliability_results[["Total"]]
    memberships <- lapply(item_names, function(v) names(factor_map)[vapply(factor_map, function(s) v %in% s, logical(1))])
    sub_item_total <- vapply(seq_along(item_names), function(i) {
      z <- memberships[[i]]
      if (length(z) != 1L) return(NA_real_)
      pos <- match(item_names[i], factor_map[[z]])
      reliability_results[[z]]$item_total[pos]
    }, numeric(1))
    sub_alpha_deleted <- vapply(seq_along(item_names), function(i) {
      z <- memberships[[i]]
      if (length(z) != 1L) return(NA_real_)
      pos <- match(item_names[i], factor_map[[z]])
      reliability_results[[z]]$alpha_deleted[pos]
    }, numeric(1))
    item_analysis <- data.frame(
      Item = item_labels,
      Variable = item_names,
      Subscale = vapply(memberships, function(z) if (length(z)) paste(z, collapse = ", ") else "", character(1)),
      `Corrected item-total r` = rr_total$item_total,
      `Alpha if deleted` = rr_total$alpha_deleted,
      `Corrected item-subscale r` = sub_item_total,
      `Subscale alpha if deleted` = sub_alpha_deleted,
      `Item mean` = vapply(X, mean, numeric(1), na.rm = TRUE),
      `Item SD` = vapply(X, stats::sd, numeric(1), na.rm = TRUE),
      check.names = FALSE,
      stringsAsFactors = FALSE
    )
  }

  # EFA ---------------------------------------------------------------------
  efa_result <- NULL
  if (isTRUE(efa)) {
    efa_result <- .r4vn_ts_efa(X, cor_method, missing, efa_method, rotation, nfactor, parallel_iter, seed)
    efa_result$loadings$Item <- item_labels[match(efa_result$loadings$Variable, item_names)]
    efa_result$loadings <- efa_result$loadings[, c("Item", "Variable", setdiff(names(efa_result$loadings), c("Item", "Variable"))), drop = FALSE]
    if (nrow(item_analysis)) {
      match_i <- match(item_analysis$Variable, efa_result$loadings$Variable)
      item_analysis$`EFA dominant factor` <- efa_result$loadings$DominantFactor[match_i]
      item_analysis$`EFA loading` <- efa_result$loadings$DominantLoading[match_i]
      item_analysis$Communality <- efa_result$loadings$Communality[match_i]
    }
  }

  # CFA ---------------------------------------------------------------------
  cfa_group_value <- NULL
  cfa_group_name <- NULL
  if (!identical(cfa_group_expr, quote(NULL))) {
    cfa_group_value <- try(eval(cfa_group_expr, envir = env), silent = TRUE)
    if (inherits(cfa_group_value, "try-error")) cfa_group_value <- NULL
    cfa_group_name <- .r4vn_ts_expr_name(cfa_group_expr, cfa_group_value, names(data), "cfa_group")
  }
  cfa_result <- NULL
  if (isTRUE(cfa)) {
    if (!is.null(cfa_group_name) && cfa_group_name %in% item_names) .r4vn_ts_stop("`cfa_group` cannot also be a scale item.")
    cfa_result <- .r4vn_ts_cfa(X, data, factor_map, ordered, estimator, cfa_missing, cfa_group_name, modification, modification_min)
    if (nrow(item_analysis)) {
      cfa_load <- cfa_result$loadings
      match_i <- match(item_analysis$Variable, cfa_load$Item)
      item_analysis$`CFA factor` <- cfa_load$Factor[match_i]
      item_analysis$`CFA loading` <- cfa_load$Loading[match_i]
      item_analysis$`CFA R2` <- cfa_load$R2[match_i]
    }
  }

  invariance_result <- NULL
  if (validity_enabled && isTRUE(invariance)) {
    if (!isTRUE(cfa)) .r4vn_ts_stop("Set `cfa = TRUE` when requesting measurement invariance.")
    invariance_result <- .r4vn_ts_invariance(X, data, factor_map, cfa_group_name,
      ordered, estimator, cfa_missing, invariance_levels)
  }

  dif_result <- data.frame()
  if (validity_enabled && isTRUE(dif)) {
    if (is.null(cfa_group_name)) .r4vn_ts_stop("DIF screening requires `cfa_group`.")
    dif_result <- .r4vn_ts_dif(X, data[[cfa_group_name]], item_labels)
  }

  # Gold-standard/criterion validity ---------------------------------------
  criterion_validity <- NULL
  gold_name <- NULL
  if (validity_enabled && !identical(gold_expr, quote(NULL))) {
    gold_value <- try(eval(gold_expr, envir = env), silent = TRUE)
    if (inherits(gold_value, "try-error")) gold_value <- NULL
    gold_name <- .r4vn_ts_expr_name(gold_expr, gold_value, names(data), "gold")
    y <- data[[gold_name]]
    nlevels_gold <- if (is.factor(y)) nlevels(droplevels(y)) else length(unique(y[!is.na(y)]))
    if (nlevels_gold == 2L) {
      rocs <- lapply(names(score_df), function(nm) .r4vn_ts_roc(y, score_df[[nm]], event, direction, conf))
      names(rocs) <- names(score_df)
      roc_table <- do.call(rbind, lapply(names(rocs), function(nm) {
        z <- rocs[[nm]]
        b <- z$best[1L, ]
        data.frame(
          Scale = nm,
          `Gold standard` = .r4vn_ts_label(data[[gold_name]], gold_name),
          `Variable name` = gold_name,
          Event = z$event,
          Direction = z$direction,
          `N event` = z$n_event,
          `N non-event` = z$n_nonevent,
          AUC = z$auc,
          `AUC p` = z$p,
          `AUC CI lower` = z$ci[1L],
          `AUC CI upper` = z$ci[2L],
          Cutoff = b$cutoff,
          Sensitivity = b$sensitivity,
          Specificity = b$specificity,
          PPV = b$ppv,
          NPV = b$npv,
          `LR+` = b$lr_pos,
          `LR-` = b$lr_neg,
          Accuracy = b$accuracy,
          DOR = b$dor,
          Youden = b$youden,
          check.names = FALSE,
          stringsAsFactors = FALSE
        )
      }))
      criterion_validity <- list(type = "roc", gold = gold_name,
        gold_label = .r4vn_ts_label(data[[gold_name]], gold_name), roc = rocs, table = roc_table)
    } else if (is.numeric(y)) {
      criterion <- do.call(rbind, lapply(names(score_df), function(nm) {
        ok <- is.finite(score_df[[nm]]) & is.finite(y)
        pearson_test <- if (sum(ok) >= 4L) try(stats::cor.test(score_df[[nm]][ok], y[ok], method = "pearson"), silent = TRUE) else NULL
        spearman_test <- if (sum(ok) >= 4L) try(suppressWarnings(stats::cor.test(score_df[[nm]][ok], y[ok], method = "spearman", exact = FALSE)), silent = TRUE) else NULL
        data.frame(
          Scale = nm,
          Criterion = .r4vn_ts_label(data[[gold_name]], gold_name),
          `Variable name` = gold_name,
          N = sum(ok),
          Pearson = if (!is.null(pearson_test) && !inherits(pearson_test, "try-error")) unname(pearson_test$estimate) else NA_real_,
          `Pearson p` = if (!is.null(pearson_test) && !inherits(pearson_test, "try-error")) pearson_test$p.value else NA_real_,
          Spearman = if (!is.null(spearman_test) && !inherits(spearman_test, "try-error")) unname(spearman_test$estimate) else NA_real_,
          `Spearman p` = if (!is.null(spearman_test) && !inherits(spearman_test, "try-error")) spearman_test$p.value else NA_real_,
          check.names = FALSE,
          stringsAsFactors = FALSE
        )
      }))
      criterion_validity <- list(type = "criterion", gold = gold_name,
        gold_label = .r4vn_ts_label(data[[gold_name]], gold_name), table = criterion)
    } else .r4vn_ts_stop("`gold` must be binary for ROC analysis or numeric for criterion correlations.")
  }

  convergent_names <- if (validity_enabled) .r4vn_ts_names(convergent, "convergent", allow_null = TRUE) else character()
  discriminant_names <- if (validity_enabled) .r4vn_ts_names(discriminant, "discriminant", allow_null = TRUE) else character()
  convergent_table <- .r4vn_ts_external_validity(score_df, data, convergent_names,
    convergent_min, cor_method, conf, "Convergent")
  discriminant_table <- .r4vn_ts_external_validity(score_df, data, discriminant_names,
    discriminant_max, cor_method, conf, "Discriminant")

  known_groups_result <- NULL
  known_groups_name <- NULL
  if (validity_enabled && !identical(known_groups_expr, quote(NULL))) {
    kg_value <- try(eval(known_groups_expr, envir = env), silent = TRUE)
    if (inherits(kg_value, "try-error")) kg_value <- NULL
    known_groups_name <- .r4vn_ts_expr_name(known_groups_expr, kg_value, names(data), "known_groups")
    known_groups_result <- .r4vn_ts_known_groups(score_df, data[[known_groups_name]],
      .r4vn_ts_label(data[[known_groups_name]], known_groups_name), conf)
  }

  # Optionally create score variables --------------------------------------
  scored_data <- data
  created <- character()
  if (!identical(name, FALSE) && !is.null(name)) {
    prefix <- if (isTRUE(name)) "scale" else as.character(name)[1L]
    if (!nzchar(prefix)) prefix <- "scale"
    for (nm in names(score_df)) {
      suffix <- if (nm == "Total") "total" else .r4vn_ts_clean_factor_name(nm)
      candidate <- make.names(paste0(prefix, "_", suffix))
      new_name <- tail(make.unique(c(names(scored_data), candidate), sep = "_"), 1L)
      scored_data[[new_name]] <- score_df[[nm]]
      attr(scored_data[[new_name]], "label") <- paste0(if (nm == "Total") "Total scale" else nm, " score (", score, ")")
      created <- c(created, new_name)
    }

    committed <- .r4vn_ts_commit_data(scored_data, data_expr, data_supplied, env)
    if (!isTRUE(committed)) .r4vn_ts_warn("Score variables were returned in `$scored_data` but could not be written back to the supplied or active data.")
  }

  # Tables for output -------------------------------------------------------
  efa_diag_table <- efa_msa_table <- efa_parallel_table <- efa_loading_table <- efa_fit_table <- efa_phi_table <- data.frame()
  if (!is.null(efa_result)) {
    d <- efa_result$diagnostics
    efa_diag_table <- data.frame(
      Statistic = c("KMO overall", "Bartlett chi-square", "Bartlett df", "Bartlett p", "Correlation determinant", "Parallel-analysis factors", "Retained factors"),
      Value = c(d$kmo, d$bartlett["chisq"], d$bartlett["df"], d$bartlett["p"], d$determinant, efa_result$parallel$suggested, efa_result$nfactor),
      stringsAsFactors = FALSE
    )
    efa_msa_table <- data.frame(Item = item_labels, Variable = item_names,
      ItemMSA = unname(d$msa[item_names]), stringsAsFactors = FALSE)
    efa_parallel_table <- efa_result$parallel$table
    efa_loading_table <- efa_result$loadings
    efa_fit_table <- data.frame(Method = efa_result$fit$method, Rotation = rotation, ChiSquare = efa_result$fit$statistic, df = efa_result$fit$df, p = efa_result$fit$p, stringsAsFactors = FALSE)
    if (!is.null(efa_result$fit$phi) && nrow(efa_result$fit$phi) > 1L) {
      phi <- efa_result$fit$phi
      efa_phi_table <- data.frame(Factor = if (is.null(rownames(phi))) paste0("F", seq_len(nrow(phi))) else rownames(phi), phi, check.names = FALSE, stringsAsFactors = FALSE)
    }
  }

  cfa_fit_table <- cfa_loading_table <- cfa_rel_table <- cfa_cor_table <- htmt_table <- mi_table <- fornell_table <- data.frame()
  if (!is.null(cfa_result)) {
    cfa_fit_table <- data.frame(Statistic = names(cfa_result$measures), Value = as.numeric(cfa_result$measures), stringsAsFactors = FALSE)
    cfa_loading_table <- cfa_result$loadings
    cfa_rel_table <- cfa_result$reliability
    cfa_cor_table <- cfa_result$factor_correlations
    htmt_table <- cfa_result$htmt
    mi_table <- cfa_result$modification
    if (nrow(cfa_rel_table) && nrow(cfa_cor_table) && !"Group" %in% names(cfa_rel_table)) {
      fornell_table <- do.call(rbind, lapply(seq_len(nrow(cfa_cor_table)), function(i) {
        f1 <- cfa_cor_table$Factor1[i]; f2 <- cfa_cor_table$Factor2[i]
        s1 <- cfa_rel_table$SqrtAVE[match(f1, cfa_rel_table$Factor)]
        s2 <- cfa_rel_table$SqrtAVE[match(f2, cfa_rel_table$Factor)]
        r <- abs(cfa_cor_table$Correlation[i])
        data.frame(Factor1=f1, Factor2=f2, `Absolute factor correlation`=r,
          `Square root AVE factor 1`=s1, `Square root AVE factor 2`=s2,
          `Fornell-Larcker supported`=ifelse(is.finite(r) & is.finite(s1) & is.finite(s2) & s1>r & s2>r,"Yes","No"),
          check.names=FALSE, stringsAsFactors=FALSE)
      }))
    }
  }

  criterion_table <- if (!is.null(criterion_validity)) criterion_validity$table else data.frame()
  reliability_score_table <- if (nrow(reliability_summary)) reliability_summary[, c("Scale","Items","ValidScores","Score mean","Score SD","Score min","Score max","Minimum valid items"),drop=FALSE] else data.frame()
  internal_consistency_table <- if (nrow(reliability_summary)) reliability_summary[, c("Scale","Alpha","Alpha CI lower","Alpha CI upper","Alpha CI method","Standardized alpha","Ordinal alpha","Omega total","Omega CI lower","Omega CI upper","Omega hierarchical","KR20","Mean inter-item r","Median inter-item r"),drop=FALSE] else data.frame()
  reliability_bounds_table <- if (nrow(reliability_summary)) reliability_summary[, c("Scale","Split-half r","Spearman-Brown","Mean split-half reliability","Best split-half reliability (lambda 4)","Guttman lambda 1","Guttman lambda 2","Guttman lambda 3","Guttman lambda 4","Guttman lambda 5","Guttman lambda 6","Greatest lower bound"),drop=FALSE] else data.frame()
  stability_table <- if (!is.null(stability)) stability$table else data.frame()
  equivalence_table <- if (!is.null(equivalence)) equivalence$table else data.frame()
  rater_table <- if (!is.null(rater_result)) rater_result$table else data.frame()
  content_item_table <- if (!is.null(content_result)) content_result$item else data.frame()
  content_summary_table <- if (!is.null(content_result)) content_result$summary else data.frame()
  known_desc_table <- if (!is.null(known_groups_result)) known_groups_result$descriptive else data.frame()
  known_test_table <- if (!is.null(known_groups_result)) known_groups_result$tests else data.frame()
  responsiveness_table <- if (!is.null(responsiveness)) responsiveness$table else data.frame()
  invariance_table <- if (!is.null(invariance_result)) invariance_result$table else data.frame()
  coverage_table <- data.frame(
    Domain = c("Internal consistency", "Stability/test-retest", "Parallel-form equivalence",
      "Inter-rater reliability", "Measurement error", "Content validity", "Face validity",
      "Structural validity (EFA)", "Structural validity (CFA)", "Convergent validity",
      "Discriminant validity", "Known-groups validity", "Criterion validity",
      "Measurement invariance", "Differential item functioning", "Responsiveness"),
    Status = c(if (isTRUE(reliability)) "Reported" else "Not requested",
      if (!is.null(stability)) "Reported" else "Not supplied",
      if (!is.null(equivalence)) "Reported" else "Not supplied",
      if (!is.null(rater_result)) "Reported" else "Not supplied",
      if (!is.null(stability) || !is.null(equivalence) || !is.null(rater_result)) "Reported" else "Requires repeated/rater data",
      if (!is.null(content_result)) "Reported" else "Not supplied",
      "Qualitative assessment; not calculated",
      if (!is.null(efa_result)) "Reported" else "Not requested",
      if (!is.null(cfa_result)) "Reported" else "Not requested",
      if (nrow(convergent_table) || nrow(cfa_rel_table)) "Reported" else "Not supplied",
      if (nrow(discriminant_table) || nrow(htmt_table) || nrow(fornell_table)) "Reported" else "Not supplied",
      if (!is.null(known_groups_result)) "Reported" else "Not supplied",
      if (!is.null(criterion_validity)) "Reported" else "Not supplied",
      if (!is.null(invariance_result)) "Reported" else "Not requested",
      if (nrow(dif_result)) "Reported" else "Not requested",
      if (!is.null(responsiveness)) "Reported" else "Not supplied"),
    stringsAsFactors = FALSE)
  tables <- list(
    `Analysis coverage` = coverage_table,
    `Item descriptive statistics` = descriptive_table,
    `Item reliability diagnostics` = item_analysis,
    `Scale score summary` = reliability_score_table,
    `Internal consistency` = internal_consistency_table,
    `Split-half and reliability lower bounds` = reliability_bounds_table,
    `Test-retest reliability and measurement error` = stability_table,
    `Parallel-form reliability` = equivalence_table,
    `Inter-rater reliability` = rater_table,
    `Content validity summary` = content_summary_table,
    `Content validity by item` = content_item_table,
    `EFA factorability` = efa_diag_table,
    `EFA item MSA` = efa_msa_table,
    `EFA model fit` = efa_fit_table,
    `EFA factor correlations` = efa_phi_table,
    `Parallel analysis` = efa_parallel_table,
    `EFA loadings` = efa_loading_table,
    `CFA fit` = cfa_fit_table,
    `CFA standardized loadings` = cfa_loading_table,
    `CFA composite reliability and AVE` = cfa_rel_table,
    `CFA factor correlations` = cfa_cor_table,
    `Fornell-Larcker discriminant validity` = fornell_table,
    `HTMT discriminant validity` = htmt_table,
    `CFA modification indices` = mi_table,
    `Measurement invariance` = invariance_table,
    `Differential item functioning screening` = dif_result,
    `Convergent validity` = convergent_table,
    `Discriminant validity` = discriminant_table,
    `Known-groups descriptive statistics` = known_desc_table,
    `Known-groups validity` = known_test_table,
    `Gold-standard or criterion validity` = criterion_table,
    Responsiveness = responsiveness_table
  )
  flat_data <- .r4vn_ts_bind_fill(tables)

  # Plot specifications are shared by the R graphics device and the Viewer.
  plot_specs <- list()
  default_plot_types <- c("correlation", "scores", "scree", "loadings",
    "cfa_loadings", "roc", "test_retest", "parallel_form", "inter_rater",
    "content", "invariance")
  add_plot <- function(name, spec) {
    if (!isTRUE(plot)) return(invisible(NULL))
    allowed <- (identical(plot_types, "auto") && name %in% default_plot_types) ||
      identical(plot_types, "all") || name %in% as.character(plot_types)
    if (allowed) plot_specs[[make.unique(c(names(plot_specs), name))[length(plot_specs)+1L]]] <<- spec
    invisible(NULL)
  }
  missing_pct <- 100 * vapply(X, function(z) mean(is.na(z)), numeric(1))
  floor_pct <- if (nrow(descriptive_table)) descriptive_table$`Floor %` else numeric()
  ceiling_pct <- if (nrow(descriptive_table)) descriptive_table$`Ceiling %` else numeric()
  add_plot("items", list(type="item_diagnostics", labels=item_labels, missing=missing_pct,
    floor=floor_pct, ceiling=ceiling_pct, title="Item completeness and floor/ceiling effects"))
  response_levels <- sort(unique(unlist(X_raw, use.names=FALSE)))
  response_levels <- response_levels[is.finite(response_levels)]
  response_prop <- sapply(X_raw, function(z) {
    tab <- table(factor(z, levels=response_levels)); if (sum(tab)) as.numeric(tab)/sum(tab) else rep(NA_real_,length(tab))
  })
  rownames(response_prop) <- paste0("Score ",response_levels)
  add_plot("distributions", list(type="item_distribution", labels=item_labels,
    proportions=response_prop, title="Item response distributions"))
  R_plot <- if (length(reliability_results)) reliability_results$Total$correlation else .r4vn_ts_cor(X, cor_method, missing)
  add_plot("correlation", list(type="correlation", matrix=R_plot, labels=item_labels, title="Item correlation matrix"))
  if (nrow(item_analysis)) add_plot("reliability", list(type="reliability", labels=item_labels,
    item_total=item_analysis$`Corrected item-total r`, alpha_deleted=item_analysis$`Alpha if deleted`,
    alpha=reliability_results$Total$alpha, title="Item reliability diagnostics"))
  add_plot("scores", list(type="scores", values=score_df, title="Scale score distributions"))
  if (!is.null(efa_result)) {
    add_plot("scree", list(type="scree", x=efa_parallel_table$Component,
      observed=efa_parallel_table$Observed, random=efa_parallel_table$Random95,
      title="Parallel-analysis scree plot"))
    L <- efa_result$fit$loadings; rownames(L) <- item_labels[match(rownames(L), item_names)]
    add_plot("loadings", list(type="loadings", loadings=L, title="Exploratory factor loadings"))
  }
  if (!is.null(cfa_result) && nrow(cfa_result$loadings)) {
    ld <- cfa_result$loadings
    ld$.factor_plot <- if ("Group" %in% names(ld)) paste(ld$Group,ld$Factor,sep=" \u2014 ") else ld$Factor
    factors_cfa <- unique(ld$.factor_plot); items_cfa <- unique(ld$`Item label`)
    Lc <- matrix(NA_real_,nrow=length(items_cfa),ncol=length(factors_cfa),dimnames=list(items_cfa,factors_cfa))
    for (i in seq_len(nrow(ld))) Lc[ld$`Item label`[i],ld$.factor_plot[i]] <- ld$Loading[i]
    add_plot("cfa_loadings",list(type="loadings",loadings=Lc,title="Confirmatory factor loadings"))
  }
  if (!is.null(criterion_validity) && criterion_validity$type == "roc") {
    curves <- lapply(criterion_validity$roc, `[[`, "curve")
    names(curves) <- paste0(names(curves), " (AUC=", .r4vn_ts_num(vapply(criterion_validity$roc, `[[`, numeric(1), "auc"), 2), ")")
    add_plot("roc", list(type="roc", curves=curves, xlim=c(0,1), ylim=c(0,1),
      title=paste0("ROC curves: ", criterion_validity$gold_label)))
  }
  if (!is.null(stability)) for (nm in names(stability$details)) {
    z <- stability$details[[nm]]; add_plot("test_retest", list(type="agreement", first=z$first, second=z$second, title=paste0("Test-retest: ", nm)))
  }
  if (!is.null(equivalence)) for (nm in names(equivalence$details)) {
    z <- equivalence$details[[nm]]; add_plot("parallel_form", list(type="agreement", first=z$first, second=z$second, title=paste0("Parallel forms: ", nm)))
  }
  if (!is.null(rater_result)) {
    vals <- c(Kappa=rater_result$table$`Unweighted/Fleiss kappa`,
      `ICC absolute`=rater_result$table$`ICC absolute single`,
      `ICC consistency`=rater_result$table$`ICC consistency single`)
    keep <- is.finite(vals); if (any(keep)) add_plot("inter_rater",
      list(type="coefficients",values=vals[keep],labels=names(vals)[keep],title="Inter-rater reliability"))
  }
  if (!is.null(content_result)) add_plot("content",
    list(type="coefficients",values=content_result$item$`I-CVI`,labels=content_result$item$Item,title="Item-level content validity"))
  external_plot <- rbind(convergent_table,discriminant_table)
  if (nrow(external_plot)) add_plot("external_validity",list(type="validity_forest",
    estimate=external_plot$Correlation,lower=external_plot$`CI lower`,upper=external_plot$`CI upper`,
    labels=paste(external_plot$Scale,external_plot$Variable,sep=" \u2014 "),title="Convergent and discriminant validity"))
  if (!is.null(known_groups_result)) for (nm in names(score_df)) add_plot("known_groups",
    list(type="known_groups", score=score_df[[nm]], group=known_groups_result$group,
      group_label=.r4vn_ts_label(data[[known_groups_name]], known_groups_name), title=paste0("Known-groups validity: ", nm)))
  if (!is.null(responsiveness)) for (nm in names(responsiveness$details)) {
    z <- responsiveness$details[[nm]]; add_plot("responsiveness",
      list(type="responsiveness", first=z$first, second=z$second, title=paste0("Responsiveness: ", nm)))
  }
  if (!is.null(invariance_result) && nrow(invariance_result$table)) add_plot("invariance",
    list(type="invariance",data=invariance_result$table,title="Measurement invariance: change in fit"))
  plot_html <- if (length(plot_specs)) paste(vapply(seq_along(plot_specs), function(i) {
    render_plot <- if (viewer_plot_format == "png") .r4vn_ts_plot_png else .r4vn_ts_plot_svg
    z <- try(render_plot(plot_specs[[i]]), silent=TRUE)
    # Keep the report usable on unusual headless systems whose PNG device is
    # unavailable; SVG remains a self-contained fallback.
    if (inherits(z, "try-error") && viewer_plot_format == "png")
      z <- try(.r4vn_ts_plot_svg(plot_specs[[i]]), silent = TRUE)
    if (inherits(z,"try-error")) "" else paste0("<div class=\"ts-chart\"><h3>",
      .r4vn_ts_escape(.r4vn_ts_or(plot_specs[[i]]$title, names(plot_specs)[i])), "</h3>", z, "</div>")
  }, character(1)), collapse="") else ""

  section_number <- 0L
  section_table <- function(df, caption) {
    if (!is.data.frame(df) || !nrow(df)) return("")
    section_number <<- section_number + 1L
    transpose <- nrow(df) <= 5L && ncol(df) >= 9L && ncol(df) > nrow(df) + 4L
    .r4vn_ts_html_table(df, paste0(section_number, ". ", caption), digit,
      p_digit, transpose = transpose)
  }

  html_sections <- c(
    section_table(coverage_table, "Analysis coverage"),
    section_table(descriptive_table, "Item descriptive statistics"),
    section_table(item_analysis, "Item reliability diagnostics"),
    section_table(reliability_score_table, "Scale and subscale score summary"),
    section_table(internal_consistency_table, "Internal consistency"),
    section_table(reliability_bounds_table, "Split-half reliability and lower-bound coefficients")
  )

  if (!is.null(efa_result)) {
    html_sections <- c(html_sections,
      section_table(efa_diag_table, "EFA factorability and retained factors"),
      section_table(efa_msa_table, "Item-level sampling adequacy"),
      section_table(efa_fit_table, "EFA model"),
      section_table(efa_phi_table, "EFA factor correlations"),
      section_table(efa_parallel_table, "Parallel analysis"),
      section_table(efa_loading_table, "Exploratory factor loadings")
    )
  }

  if (!is.null(cfa_result)) {
    html_sections <- c(html_sections,
      section_table(cfa_fit_table, "Confirmatory factor analysis fit"),
      paste0("<div class=\"ts-model\"><strong>CFA model</strong><pre>", .r4vn_ts_escape(cfa_result$model), "</pre><div>Estimator: ", .r4vn_ts_escape(cfa_result$estimator), "</div></div>"),
      section_table(cfa_loading_table, "CFA standardized loadings"),
      section_table(cfa_rel_table, "Composite reliability and convergent validity"),
      section_table(cfa_cor_table, "Factor correlations"),
      section_table(fornell_table, "Fornell-Larcker discriminant validity"),
      section_table(htmt_table, "HTMT discriminant validity"),
      section_table(mi_table, "Modification indices")
    )
  }

  html_sections <- c(html_sections,
    section_table(stability_table, "Test-retest reliability, agreement, and measurement error"),
    section_table(equivalence_table, "Parallel-form reliability and agreement"),
    section_table(rater_table, "Inter-rater reliability and agreement"),
    section_table(content_summary_table, "Content validity summary"),
    section_table(content_item_table, "Content validity by item"),
    section_table(invariance_table, "Measurement invariance"),
    section_table(dif_result, "Differential item functioning screening"),
    section_table(convergent_table, "Convergent validity"),
    section_table(discriminant_table, "Discriminant validity"),
    section_table(known_desc_table, "Known-groups descriptive statistics"),
    section_table(known_test_table, "Known-groups validity"),
    section_table(criterion_table, if (!is.null(criterion_validity) && criterion_validity$type == "roc") "Gold-standard ROC analysis" else "Criterion validity"),
    section_table(responsiveness_table, "Responsiveness"),
    plot_html)

  interpretation_text <- character()
  if (isTRUE(interpretation)) {
    alpha_total <- if (nrow(reliability_summary)) reliability_summary$Alpha[reliability_summary$Scale == "Total"][1L] else NA_real_
    if (is.finite(alpha_total)) interpretation_text <- c(interpretation_text,
      sprintf("The total-scale Cronbach alpha was %.2f; interpret this together with dimensionality, item content, and the intended use of the score.", alpha_total))
    if (!is.null(efa_result)) interpretation_text <- c(interpretation_text,
      sprintf("The overall KMO was %.2f, Bartlett's test p was %s, and parallel analysis suggested %d factor(s).",
        efa_result$diagnostics$kmo, .r4vn_ts_p(efa_result$diagnostics$bartlett["p"], p_digit), efa_result$parallel$suggested))
    if (!is.null(cfa_result)) interpretation_text <- c(interpretation_text,
      sprintf("CFA fit: CFI %.3f, TLI %.3f, RMSEA %.3f, and SRMR %.3f; prespecified criteria and substantive plausibility should guide conclusions.",
        cfa_result$measures["cfi"], cfa_result$measures["tli"], cfa_result$measures["rmsea"], cfa_result$measures["srmr"]))
    if (nrow(known_test_table)) interpretation_text <- c(interpretation_text,
      paste0("Known-groups comparisons were produced for ", nrow(known_test_table), " scale score(s); assess whether effect direction and magnitude match prespecified hypotheses."))
    if (!is.null(criterion_validity) && criterion_validity$type == "roc") interpretation_text <- c(interpretation_text,
      "ROC estimates describe discrimination in this sample and should be externally validated before a cutoff is used clinically.")
  }
  if (length(interpretation_text)) html_sections <- c(html_sections,
    paste0("<div class=\"ts-interpretation\"><h3>Interpretation</h3>",
      paste0("<p>", .r4vn_ts_escape(interpretation_text), "</p>", collapse=""), "</div>"))

  negative_items <- if (nrow(item_analysis)) item_analysis$Variable[is.finite(item_analysis$`Corrected item-total r`) & item_analysis$`Corrected item-total r` < 0] else character()
  notes <- c(
    "Categorical responses are shown as n (%). Item means, reliability, EFA, CFA, subscale scores, and total scores use reverse-scored values where specified.",
    "Corrected item-total correlation excludes the item from its comparison total. Omega total is estimated from a one-factor principal-axis solution.",
    if (length(reverse_names)) paste0("Reverse-scored items: ", paste(reverse_names, collapse = ", "), ".") else NULL,
    if (length(negative_items)) paste0("Items with negative corrected item-total correlation: ", paste(negative_items, collapse = ", "), ". Check coding and construct direction.") else NULL,
    if (length(created)) paste0("Created score variables: ", paste(created, collapse = ", "), ".") else NULL,
    if (!is.null(criterion_validity) && criterion_validity$type == "roc") "AUC confidence intervals use the Hanley-McNeil large-sample approximation; the displayed cutoff maximizes Youden's index." else NULL,
    "Face validity is a qualitative judgement by target users and experts; tabscale records it in the coverage table but does not manufacture a numeric coefficient.",
    "Ordinal alpha, omega hierarchical, and the greatest lower bound are reported when the optional psych package is available; blank values are not interpreted as zero.",
    if (nrow(dif_result)) "DIF results are screening analyses based on nested linear models and should be confirmed with an ordinal-logistic or item-response model when items are ordinal." else NULL,
    if (isTRUE(cfa) && length(cfa_result$ordered)) "Ordinal CFA was estimated with the declared ordered items." else NULL
  )
  note_html <- paste0("<div class=\"ts-notes\">", paste0("<div>", .r4vn_ts_escape(notes), "</div>", collapse = ""), "</div>")

  title <- .r4vn_ts_or(title, "Scale analysis")
  title_html <- paste0("<h2>", .r4vn_ts_escape(title), "</h2>")
  table_block <- paste0("<section class=\"r4vn-table r4vn-tabscale\">", title_html, paste(html_sections[nzchar(html_sections)], collapse = ""), note_html, "</section>")

  css <- switch(template,
    journal = "body{font-family:Arial,sans-serif;color:#111;margin:20px}.r4vn-tabscale{max-width:100%}h2{font-size:20px;margin:0 0 14px}h3{font-size:15px;margin:20px 0 7px}.ts-section{max-width:100%;overflow-x:auto}table{border-collapse:collapse;font-size:12px;min-width:620px}th{border-top:2px solid #222;border-bottom:1px solid #555;padding:6px 8px;text-align:right;white-space:nowrap}th:first-child,td:first-child{text-align:left}td{border-bottom:1px solid #ddd;padding:5px 8px;text-align:right;vertical-align:top;white-space:nowrap}.ts-notes{font-size:11px;line-height:1.4;margin-top:14px}.ts-chart{max-width:900px;overflow:hidden;margin:18px 0}.ts-chart img.ts-plot-image,.ts-chart svg{display:block;width:100%;height:auto;max-width:900px}.ts-chart svg text{font-family:Arial,sans-serif!important;letter-spacing:normal!important;word-spacing:normal!important}.ts-interpretation{border-left:4px solid #555;padding:2px 12px;margin:18px 0;font-size:12px}.ts-svg{width:100%;height:auto;color:#222;background:#fff}.ts-model{font-size:12px;margin:12px 0}.ts-model pre{background:#f6f6f6;padding:8px;white-space:pre-wrap}.table-separator{height:26px}",
    clean = "body{font-family:system-ui,Arial,sans-serif;color:#222;margin:20px}.r4vn-tabscale{max-width:100%}h2{font-size:21px}h3{font-size:15px;margin-top:20px}.ts-section{overflow-x:auto}table{border-collapse:separate;border-spacing:0;font-size:12px;min-width:620px;border:1px solid #ddd;border-radius:5px}th{background:#f5f5f5;font-weight:600}th,td{padding:6px 8px;border-bottom:1px solid #e4e4e4;text-align:right;white-space:nowrap}th:first-child,td:first-child{text-align:left}.ts-notes{font-size:11px;line-height:1.45;margin-top:14px}.ts-chart{max-width:900px;overflow:hidden;margin:18px 0}.ts-chart img.ts-plot-image,.ts-chart svg{display:block;width:100%;height:auto;max-width:900px}.ts-chart svg text{font-family:Arial,sans-serif!important;letter-spacing:normal!important;word-spacing:normal!important}.ts-svg{width:100%;max-width:900px;color:#333}.ts-model pre{background:#f5f5f5;padding:8px}.table-separator{height:26px}",
    minimal = "body{font-family:Arial,sans-serif;margin:16px}h2{font-size:18px}h3{font-size:14px}.ts-section{overflow-x:auto}table{border-collapse:collapse;font-size:11px}th,td{padding:4px 6px;border-bottom:1px solid #ccc;text-align:right;white-space:nowrap}th:first-child,td:first-child{text-align:left}.ts-notes{font-size:10px;margin-top:10px}.ts-chart{max-width:900px;overflow:hidden;margin:16px 0}.ts-chart img.ts-plot-image,.ts-chart svg{display:block;width:100%;height:auto;max-width:900px}.ts-chart svg text{font-family:Arial,sans-serif!important;letter-spacing:normal!important;word-spacing:normal!important}.ts-svg{width:100%;max-width:900px;color:#111}.table-separator{height:20px}"
  )

  blocks <- table_block
  if (inherits(append, "r4vn_tab") && !is.null(append$blocks)) blocks <- c(append$blocks, table_block)
  document <- paste0("<!DOCTYPE html><html><head><meta charset=\"UTF-8\"><meta name=\"viewport\" content=\"width=device-width,initial-scale=1\"><style>", css, "</style></head><body>", paste(blocks, collapse = "<div class=\"table-separator\"></div>"), "</body></html>")

  if (is.null(file)) file <- tempfile(pattern = "r4vn-tabscale-", fileext = ".html")
  if (!is.character(file) || length(file) != 1L || !nzchar(file)) .r4vn_ts_stop("`file` must be one valid path.")
  if (!grepl("[.]html?$", file, ignore.case = TRUE)) file <- paste0(file, ".html")
  if (is.character(append) && length(append) == 1L && file.exists(append)) {
    old <- paste(readLines(append, warn = FALSE, encoding = "UTF-8"), collapse = "\n")
    if (grepl("</body>", old, fixed = TRUE)) {
      document <- sub("</body>", paste0("<div class=\"table-separator\"></div>", table_block, "</body>"), old, fixed = TRUE)
      file <- append
    }
  }
  writeLines(enc2utf8(document), file, useBytes = TRUE)

  output <- list(
    report = report,
    data = flat_data,
    tables = tables,
    coverage = coverage_table,
    descriptive = descriptive_table,
    item_analysis = item_analysis,
    reliability = if (isTRUE(raw)) reliability_results else reliability_summary,
    reliability_summary = reliability_summary,
    scores = score_df,
    score_details = score_objects,
    scored_data = scored_data,
    created = created,
    items = item_names,
    reverse = reverse_names,
    factor = factor_map,
    efa = efa_result,
    cfa = cfa_result,
    stability = stability,
    test_retest = stability,
    parallel_form = equivalence,
    inter_rater = rater_result,
    content_validity = content_result,
    convergent_validity = convergent_table,
    discriminant_validity = discriminant_table,
    known_groups = known_groups_result,
    invariance = invariance_result,
    dif = dif_result,
    responsiveness = responsiveness,
    validity = criterion_validity,
    criterion_validity = criterion_validity,
    validity_all = list(content=content_result, structural=list(efa=efa_result,cfa=cfa_result),
      convergent=convergent_table, discriminant=list(external=discriminant_table,htmt=htmt_table,fornell_larcker=fornell_table),
      known_groups=known_groups_result, criterion=criterion_validity, invariance=invariance_result,
      dif=dif_result, responsiveness=responsiveness),
    gold = gold_name,
    plots = plot_specs,
    interpretation = interpretation_text,
    html = document,
    table_html = table_block,
    blocks = blocks,
    file = normalizePath(file, winslash = "/", mustWork = TRUE),
    call = match.call()
  )
  class(output) <- c("r4vn_tabscale", "r4vn_tab")
  if (isTRUE(plot) && length(plot_specs) && interactive()) graphics::plot(output)
  if (isTRUE(show) && interactive()) {
    viewer <- getOption("viewer")
    if (is.function(viewer)) viewer(output$file) else utils::browseURL(output$file)
  }
  invisible(output)
}

#' Print or reopen a scale analysis
#'
#' @param x An object created by `tabscale()`.
#' @param ... Additional arguments currently ignored.
#' @return The input object, invisibly.
#' @method print r4vn_tabscale
#' @export
print.r4vn_tabscale <- function(x, ...) {
  if (!inherits(x, "r4vn_tabscale")) .r4vn_ts_stop("`x` must be created by tabscale().")
  viewer <- getOption("viewer")
  if (is.function(viewer)) viewer(x$file) else utils::browseURL(x$file)
  invisible(x)
}

#' Plot diagnostics from a scale analysis
#'
#' Displays every available tabscale graphic in the R graphics device. In
#' RStudio, use the back/forward arrows in the Plots pane to review them.
#'
#' @param x An object created by `tabscale()`.
#' @param which Optional plot names or indices. The default displays all plots.
#' @param ... Additional arguments currently ignored.
#' @return The input object, invisibly.
#' @method plot r4vn_tabscale
#' @export
plot.r4vn_tabscale <- function(x, which = NULL, ...) {
  if (!inherits(x, "r4vn_tabscale")) .r4vn_ts_stop("`x` must be created by tabscale().")
  specs <- x$plots
  if (!length(specs)) {
    .r4vn_ts_warn("No plots are stored in this result. Run tabscale(..., plot = TRUE).")
    return(invisible(x))
  }
  if (!is.null(which)) {
    if (is.numeric(which)) specs <- specs[intersect(as.integer(which), seq_along(specs))]
    else specs <- specs[names(specs) %in% as.character(which) | sub("[.][0-9]+$", "", names(specs)) %in% as.character(which)]
  }
  for (spec in specs) .r4vn_ts_plot_draw(spec)
  invisible(x)
}

attr(tabscale, "r4vn_version") <- "tabscale-2.0.2-2026-09-02"