R/tablearn.R

In R4VN: Health Data Analysis and Publication-Ready Reporting

# tablearn.R
# Comprehensive learning-curve analysis for R4VN
# Designed to work with minimal hard dependencies. Standard analysis, Viewer
# reporting, and base-R plots require no add-on package. ggplot2 is used only
# when it is already installed to provide the advanced plotting path.

# Avoid R CMD check notes for ggplot2 NSE variables used internally.
if (getRversion() >= "2.15.1") {
  utils::globalVariables(c(
    ".case", ".order", ".operator", ".value", ".lower", ".upper",
    ".series", ".phase", ".current", ".window_label", ".x", ".n", ".bp",
    ".label", ".xmin", ".xmax", ".xmid", ".phase_fill", ".limit", ".signal", ".proficient", ".prof", "estimate", "phase"
  ))
}

.tablearn_or <- function(x, y) if (is.null(x)) y else x

.tablearn_stop <- function(...) stop(..., call. = FALSE)
.tablearn_warn <- function(...) warning(..., call. = FALSE)

.tablearn_lang <- function(lang, en, vi) {
  if (identical(lang, "vi")) vi else en
}

.tablearn_active_data <- function(data = NULL) {
  if (!is.null(data)) {
    if (!is.data.frame(data)) .tablearn_stop("`data` must be a data.frame.")
    return(data)
  }

  # Compatibility hooks for R4VN active-data implementations.
  f <- get0(".r4vn_get_active_data", mode = "function", inherits = TRUE)
  if (!is.null(f)) {
    ans <- try(f(), silent = TRUE)
    if (!inherits(ans, "try-error") && is.data.frame(ans)) return(ans)
  }

  opt <- getOption("R4VN.active_data", NULL)
  if (is.data.frame(opt)) return(opt)
  if (is.character(opt) && length(opt) == 1L && exists(opt, envir = .GlobalEnv, inherits = FALSE)) {
    ans <- get(opt, envir = .GlobalEnv, inherits = FALSE)
    if (is.data.frame(ans)) return(ans)
  }

  for (nm in c(".r4vn_active_data", ".R4VN_active_data", "R4VN.active.data")) {
    if (exists(nm, envir = .GlobalEnv, inherits = FALSE)) {
      ans <- get(nm, envir = .GlobalEnv, inherits = FALSE)
      if (is.data.frame(ans)) return(ans)
    }
  }

  .tablearn_stop(
    "No `data` was supplied and no R4VN active data frame could be found. ",
    "Supply `data = ...` or activate a data frame first."
  )
}

.tablearn_extract_vars_from_expr <- function(expr, data, env, allow_null = TRUE) {
  if (is.null(expr) || identical(expr, quote(NULL))) {
    if (allow_null) return(NULL)
    .tablearn_stop("A variable is required.")
  }

  # Bare variable.
  if (is.symbol(expr)) {
    nm <- as.character(expr)
    if (nm %in% names(data)) return(nm)
    val <- try(eval(expr, envir = env), silent = TRUE)
    if (!inherits(val, "try-error") && is.character(val)) return(as.character(val))
    .tablearn_stop("Variable not found in `data`: ", nm, ".")
  }

  # Character vector or an object that evaluates to character names.
  val <- try(eval(expr, envir = env), silent = TRUE)
  if (!inherits(val, "try-error")) {
    if (is.character(val)) return(as.character(val))
    if (is.numeric(val) && all(val %in% seq_along(data))) return(names(data)[val])
    # Native R4VN vars(...) objects. Resolve deferred selectors/prefixes against
    # the actual data frame when possible; this makes tablearn() consistent with
    # the rest of R4VN (e.g. vars(c.age, sex), vars(`score*`)).
    if (inherits(val, "r4vn_vars") || (is.data.frame(val) && "variable" %in% names(val))) {
      resolver <- get0(".r4vn_resolve_vars", mode = "function", inherits = TRUE)
      if (is.function(resolver)) {
        resolved <- try(resolver(val, data, default_type = "auto", strict = TRUE), silent = TRUE)
        if (!inherits(resolved, "try-error") && is.data.frame(resolved) && "variable" %in% names(resolved)) {
          return(unique(as.character(resolved$variable)))
        }
      }
      return(unique(as.character(val$variable)))
    }
  }

  # R4VN vars(...) syntax: capture the arguments directly if evaluation is unavailable.
  if (is.call(expr) && identical(as.character(expr[[1L]]), "vars")) {
    args <- as.list(expr)[-1L]
    out <- character()
    for (a in args) {
      if (is.symbol(a)) {
        out <- c(out, as.character(a))
      } else if (is.character(a)) {
        out <- c(out, a)
      } else {
        av <- try(eval(a, envir = env), silent = TRUE)
        if (!inherits(av, "try-error") && is.character(av)) out <- c(out, av)
      }
    }
    if (length(out)) return(out)
  }

  .tablearn_stop("Could not resolve variable specification: ", paste(deparse(expr), collapse = " "))
}

.tablearn_check_vars <- function(vars, data, arg = "variable") {
  if (is.null(vars)) return(invisible(TRUE))
  miss <- setdiff(vars, names(data))
  if (length(miss)) .tablearn_stop(
    "Variables not found in `data` for `", arg, "`: ", paste(miss, collapse = ", "), "."
  )
  invisible(TRUE)
}

.tablearn_detect_type <- function(y, type = "auto", event = NULL) {
  if (!identical(type, "auto")) return(type)
  yy <- y[!is.na(y)]
  if (!length(yy)) .tablearn_stop("The outcome contains no non-missing observations.")

  # Non-numeric two-level variables are unambiguously categorical.
  if (is.logical(yy) || is.factor(yy) || is.character(yy)) {
    if (length(unique(yy)) == 2L) return("binary")
  }

  if (is.numeric(yy) || is.integer(yy)) {
    u <- unique(as.numeric(yy))

    # Numeric outcomes need a conservative rule. A genuinely continuous
    # measurement can have only two observed values in a small or deliberately
    # constructed dataset (e.g. 50 and 100 minutes). Treating every two-valued
    # numeric variable as binary would silently turn means into proportions and
    # can produce false target/proficiency signals.
    #
    # Auto-detect only canonical 0/1 coding as binary. Other two-level numeric
    # codings are binary when the user explicitly identifies the event level;
    # otherwise they remain continuous. Users can always force the intended
    # interpretation with `type = "binary"` or `type = "continuous"`.
    if (length(u) == 2L) {
      if (all(sort(u) == c(0, 1))) return("binary")
      if (!is.null(event) && length(event) == 1L && !is.na(event)) {
        event_num <- suppressWarnings(as.numeric(event))
        if (is.finite(event_num) && any(u == event_num)) return("binary")
      }
      return("continuous")
    }
    return("continuous")
  }
  .tablearn_stop("Unsupported outcome type. Use numeric, integer, logical, factor, or character outcome.")
}

.tablearn_binary <- function(y, event = NULL) {
  yy <- y[!is.na(y)]
  lev <- unique(yy)
  if (length(lev) != 2L) .tablearn_stop("Binary outcome must have exactly two non-missing levels.")

  if (is.null(event)) {
    if (is.factor(y)) {
      lv <- levels(y)
      lv <- lv[lv %in% as.character(yy)]
      event <- if (length(lv) >= 2L) lv[2L] else as.character(lev[2L])
    } else if (is.logical(y)) {
      event <- TRUE
    } else if (is.numeric(y) || is.integer(y)) {
      event <- max(lev)
    } else {
      event <- as.character(lev[2L])
    }
  }

  z <- if (is.factor(y) || is.character(y)) {
    as.integer(as.character(y) == as.character(event))
  } else {
    as.integer(y == event)
  }
  z[is.na(y)] <- NA_integer_
  list(value = z, event = event)
}

.tablearn_weight_vector <- function(n, method = "equal", decay = 0.85) {
  if (n <= 0L) return(numeric())
  if (method == "equal") return(rep(1, n))
  if (method == "linear") return(seq_len(n))
  if (method == "exponential") {
    if (!is.numeric(decay) || length(decay) != 1L || decay <= 0 || decay > 1) {
      .tablearn_stop("`window_decay` must be in (0, 1] for exponential weighting.")
    }
    return(decay ^ rev(seq.int(0, n - 1L)))
  }
  rep(1, n)
}

.tablearn_weighted_mean <- function(x, w) {
  ok <- is.finite(x) & is.finite(w)
  if (!any(ok)) return(NA_real_)
  sum(x[ok] * w[ok]) / sum(w[ok])
}

.tablearn_weighted_sd <- function(x, w) {
  ok <- is.finite(x) & is.finite(w)
  x <- x[ok]; w <- w[ok]
  if (length(x) < 2L || sum(w) <= 0) return(NA_real_)
  m <- sum(w * x) / sum(w)
  # Reliability-weight approximation.
  denom <- sum(w) - sum(w^2) / sum(w)
  if (denom <= 0) return(NA_real_)
  sqrt(sum(w * (x - m)^2) / denom)
}

.tablearn_effective_n <- function(w) {
  if (!length(w) || sum(w^2) == 0) return(NA_real_)
  sum(w)^2 / sum(w^2)
}

.tablearn_wilson <- function(x, n, level = 0.95) {
  if (!is.finite(n) || n <= 0) return(c(NA_real_, NA_real_))
  p <- x / n
  z <- stats::qnorm(1 - (1 - level) / 2)
  den <- 1 + z^2 / n
  ctr <- (p + z^2 / (2 * n)) / den
  half <- z * sqrt(p * (1 - p) / n + z^2 / (4 * n^2)) / den
  c(max(0, ctr - half), min(1, ctr + half))
}

.tablearn_window_ci <- function(x, estimate, type, method, level, weights = NULL, fun = "mean") {
  x <- x[!is.na(x)]
  n <- length(x)
  if (!n) return(c(lower = NA_real_, upper = NA_real_))
  if (is.null(weights)) weights <- rep(1, n)
  weights <- weights[seq_len(min(length(weights), n))]

  if (type == "binary") {
    if (!identical(weights, rep(1, length(weights))) && method %in% c("auto", "wilson", "exact")) {
      method <- "normal"
    }
    if (method == "auto") method <- "wilson"
    if (method == "wilson") {
      ci <- .tablearn_wilson(sum(x), n, level)
      return(c(lower = ci[1], upper = ci[2]))
    }
    if (method == "exact") {
      bt <- stats::binom.test(sum(x), n, conf.level = level)
      return(c(lower = unname(bt$conf.int[1]), upper = unname(bt$conf.int[2])))
    }
    ne <- .tablearn_effective_n(weights)
    z <- stats::qnorm(1 - (1 - level) / 2)
    se <- sqrt(max(estimate * (1 - estimate), 0) / ne)
    return(c(lower = max(0, estimate - z * se), upper = min(1, estimate + z * se)))
  }

  if (fun == "median") {
    # Distribution-free sign-test style interval for the population median.
    xs <- sort(x)
    a <- (1 - level) / 2
    k <- stats::qbinom(a, n, 0.5)
    lo_i <- max(1L, as.integer(k + 1L))
    hi_i <- min(n, as.integer(n - k))
    return(c(lower = xs[lo_i], upper = xs[hi_i]))
  }

  if (method == "auto") method <- if (n < 30L) "t" else "normal"
  sdv <- .tablearn_weighted_sd(x, weights)
  ne <- .tablearn_effective_n(weights)
  se <- sdv / sqrt(ne)
  if (!is.finite(se)) return(c(lower = NA_real_, upper = NA_real_))
  q <- if (method == "t") stats::qt(1 - (1 - level) / 2, df = max(1, floor(ne - 1))) else stats::qnorm(1 - (1 - level) / 2)
  c(lower = estimate - q * se, upper = estimate + q * se)
}

.tablearn_window_index <- function(n, size, type, align, step, complete) {
  if (n <= 0L) return(data.frame(start = integer(), end = integer()))
  if (type == "none") return(data.frame(start = seq_len(n), end = seq_len(n)))
  if (type == "cumulative") {
    ends <- seq.int(1L, n, by = step)
    if (tail(ends, 1L) != n) ends <- c(ends, n)
    return(data.frame(start = rep(1L, length(ends)), end = ends))
  }
  if (type == "block") {
    starts <- seq.int(1L, n, by = step)
    ends <- pmin(starts + size - 1L, n)
    keep <- if (complete) (ends - starts + 1L) == size else rep(TRUE, length(starts))
    return(data.frame(start = starts[keep], end = ends[keep]))
  }

  # rolling
  if (align == "right") {
    ends <- seq.int(if (complete) size else 1L, n, by = step)
    starts <- pmax(1L, ends - size + 1L)
  } else if (align == "left") {
    starts <- seq.int(1L, if (complete) max(1L, n - size + 1L) else n, by = step)
    ends <- pmin(n, starts + size - 1L)
  } else {
    centers <- seq.int(1L, n, by = step)
    left <- floor((size - 1L) / 2L)
    starts <- centers - left
    ends <- starts + size - 1L
    if (complete) {
      keep <- starts >= 1L & ends <= n
      starts <- starts[keep]; ends <- ends[keep]
    } else {
      starts <- pmax(1L, starts); ends <- pmin(n, ends)
    }
  }
  data.frame(start = as.integer(starts), end = as.integer(ends))
}

.tablearn_window_fun <- function(x, outcome_type, fun, trim, weights) {
  ok <- !is.na(x)
  x <- x[ok]
  weights <- weights[ok]
  if (!length(x)) return(NA_real_)

  if (fun == "auto") fun <- if (outcome_type == "binary") "proportion" else "mean"
  if (fun %in% c("mean", "proportion", "rate")) return(.tablearn_weighted_mean(as.numeric(x), weights))
  if (fun == "median") return(stats::median(x, na.rm = TRUE))
  if (fun == "trimmed_mean") return(mean(x, trim = trim, na.rm = TRUE))
  .tablearn_stop("Unsupported `window_fun`: ", fun)
}

.tablearn_compute_windows <- function(y, x_order, outcome_type, size = 5L,
                                      type = "rolling", align = "right", step = 1L,
                                      complete = TRUE, fun = "auto", trim = 0.10,
                                      weight = "equal", decay = 0.85,
                                      ci = TRUE, ci_level = 0.95,
                                      ci_method = "auto", interval = "ci") {
  n <- length(y)
  idx <- .tablearn_window_index(n, size, type, align, step, complete)
  if (!nrow(idx)) return(data.frame())

  rows <- vector("list", nrow(idx))
  for (i in seq_len(nrow(idx))) {
    ii <- idx$start[i]:idx$end[i]
    yy <- y[ii]
    w <- .tablearn_weight_vector(length(ii), weight, decay)
    est <- .tablearn_window_fun(yy, outcome_type, fun, trim, w)
    lower <- upper <- NA_real_

    resolved_fun <- if (fun == "auto") if (outcome_type == "binary") "proportion" else "mean" else fun
    if (interval == "ci" && isTRUE(ci)) {
      ci0 <- .tablearn_window_ci(yy, est, outcome_type, ci_method, ci_level, w, resolved_fun)
      lower <- ci0[1]; upper <- ci0[2]
    } else if (interval == "sd" && outcome_type != "binary") {
      s <- .tablearn_weighted_sd(yy, w)
      lower <- est - s; upper <- est + s
    } else if (interval == "iqr" && outcome_type != "binary") {
      q <- stats::quantile(yy, c(.25, .75), na.rm = TRUE, names = FALSE, type = 7)
      lower <- q[1]; upper <- q[2]
    } else if (interval == "range" && outcome_type != "binary") {
      lower <- min(yy, na.rm = TRUE); upper <- max(yy, na.rm = TRUE)
    }

    if (align == "right") {
      xpos <- idx$end[i]
    } else if (align == "left") {
      xpos <- idx$start[i]
    } else {
      xpos <- (idx$start[i] + idx$end[i]) / 2
    }

    xo <- x_order[ii]
    xord <- if (inherits(x_order, "Date")) {
      if (align == "right") max(xo, na.rm = TRUE) else if (align == "left") min(xo, na.rm = TRUE) else as.Date(mean(as.numeric(xo), na.rm = TRUE), origin = "1970-01-01")
    } else if (is.numeric(x_order)) {
      if (align == "right") xo[length(xo)] else if (align == "left") xo[1L] else mean(xo, na.rm = TRUE)
    } else {
      xo[length(xo)]
    }

    rows[[i]] <- data.frame(
      .window_id = i,
      .start = idx$start[i],
      .end = idx$end[i],
      .x_case = xpos,
      .n = sum(!is.na(yy)),
      .events = if (outcome_type == "binary") sum(yy, na.rm = TRUE) else NA_real_,
      .value = est,
      .lower = lower,
      .upper = upper,
      .window_label = paste0(idx$start[i], "-", idx$end[i]),
      stringsAsFactors = FALSE
    )
    rows[[i]]$.x_order <- xord
  }
  out <- do.call(rbind, rows)
  rownames(out) <- NULL
  out
}

.tablearn_ewma <- function(y, lambda = 0.2, init = "first", target = NULL) {
  if (!is.numeric(lambda) || length(lambda) != 1L || lambda <= 0 || lambda > 1) {
    .tablearn_stop("`ewma_lambda` must be in (0, 1].")
  }
  n <- length(y)
  if (!n) return(data.frame())
  z <- rep(NA_real_, n)
  ynum <- as.numeric(y)
  first_ok <- which(!is.na(ynum))[1L]
  if (is.na(first_ok)) return(data.frame(.case = seq_len(n), .value = z))

  z0 <- switch(init,
    first = ynum[first_ok],
    mean = mean(ynum, na.rm = TRUE),
    target = {
      if (is.null(target) || !is.numeric(target) || length(target) != 1L) .tablearn_stop("`target` is required when `ewma_init = 'target'`.")
      target
    },
    ynum[first_ok]
  )
  prev <- z0
  for (i in seq_len(n)) {
    if (is.na(ynum[i])) {
      z[i] <- prev
    } else {
      prev <- lambda * ynum[i] + (1 - lambda) * prev
      z[i] <- prev
    }
  }
  data.frame(.case = seq_len(n), .value = z)
}

.tablearn_smooth <- function(x, y, method = "loess", span = 0.6, df = NULL,
                             ci = TRUE, ci_level = 0.95, outcome_type = "continuous") {
  ok <- is.finite(x) & is.finite(y)
  x <- x[ok]; y <- y[ok]
  if (length(x) < 4L || method == "none") return(data.frame())
  ord <- order(x); x <- x[ord]; y <- y[ord]
  gx <- sort(unique(x))
  z <- stats::qnorm(1 - (1 - ci_level) / 2)

  if (method == "loess") {
    fit <- try(stats::loess(y ~ x, span = span, na.action = stats::na.exclude, control = stats::loess.control(surface = "direct")), silent = TRUE)
    if (inherits(fit, "try-error")) return(data.frame())
    pr <- try(stats::predict(fit, newdata = data.frame(x = gx), se = isTRUE(ci)), silent = TRUE)
    if (inherits(pr, "try-error")) return(data.frame())
    if (is.list(pr)) {
      val <- as.numeric(pr$fit); se <- as.numeric(pr$se.fit)
    } else {
      val <- as.numeric(pr); se <- rep(NA_real_, length(val))
    }
    lo <- val - z * se; hi <- val + z * se
    if (outcome_type == "binary") {
      val <- pmin(1, pmax(0, val)); lo <- pmin(1, pmax(0, lo)); hi <- pmin(1, pmax(0, hi))
    }
    return(data.frame(.case = gx, .value = val, .lower = lo, .upper = hi))
  }

  if (method == "lm") {
    fit <- if (outcome_type == "binary") stats::glm(y ~ x, family = stats::binomial()) else stats::lm(y ~ x)
    pr <- if (outcome_type == "binary") {
      p <- stats::predict(fit, newdata = data.frame(x = gx), type = "link", se.fit = TRUE)
      val <- stats::plogis(p$fit)
      lo <- stats::plogis(p$fit - z * p$se.fit); hi <- stats::plogis(p$fit + z * p$se.fit)
      list(val = val, lo = lo, hi = hi)
    } else {
      p <- stats::predict(fit, newdata = data.frame(x = gx), se.fit = TRUE)
      list(val = as.numeric(p$fit), lo = as.numeric(p$fit - z * p$se.fit), hi = as.numeric(p$fit + z * p$se.fit))
    }
    return(data.frame(.case = gx, .value = pr$val, .lower = pr$lo, .upper = pr$hi))
  }

  if (method == "spline") {
    fit <- try(stats::smooth.spline(x, y, df = df), silent = TRUE)
    if (inherits(fit, "try-error")) return(data.frame())
    pr <- stats::predict(fit, gx)
    vv <- as.numeric(pr$y)
    if (outcome_type == "binary") vv <- pmin(1, pmax(0, vv))
    return(data.frame(.case = gx, .value = vv, .lower = NA_real_, .upper = NA_real_))
  }

  if (method == "gam") {
    if (!requireNamespace("mgcv", quietly = TRUE)) {
      .tablearn_warn("`smooth_method = 'gam'` requires package `mgcv`; falling back to loess.")
      return(.tablearn_smooth(x, y, "loess", span, df, ci, ci_level, outcome_type))
    }
    k <- max(3L, min(10L, length(unique(x)) - 1L))
    s <- mgcv::s
    fit <- if (outcome_type == "binary") {
      mgcv::gam(y ~ s(x, k = k), family = stats::binomial())
    } else {
      mgcv::gam(y ~ s(x, k = k))
    }
    p <- stats::predict(fit, newdata = data.frame(x = gx), type = "link", se.fit = TRUE)
    if (outcome_type == "binary") {
      val <- stats::plogis(p$fit); lo <- stats::plogis(p$fit - z * p$se.fit); hi <- stats::plogis(p$fit + z * p$se.fit)
    } else {
      val <- as.numeric(p$fit); lo <- val - z * p$se.fit; hi <- val + z * p$se.fit
    }
    return(data.frame(.case = gx, .value = val, .lower = lo, .upper = hi))
  }

  data.frame()
}

.tablearn_fit_hinge <- function(x, y, bp, outcome_type) {
  hinge <- pmax(0, x - bp)
  dat <- data.frame(y = y, x = x, hinge = hinge)
  fit <- try(
    if (outcome_type == "binary") stats::glm(y ~ x + hinge, family = stats::binomial(), data = dat)
    else if (outcome_type == "count") stats::glm(y ~ x + hinge, family = stats::poisson(), data = dat)
    else stats::lm(y ~ x + hinge, data = dat),
    silent = TRUE
  )
  if (inherits(fit, "try-error")) return(NULL)
  fit
}

.tablearn_breakpoint <- function(x, y, outcome_type, min_n = 8L, grid = NULL, boot = 0L, ci_level = 0.95) {
  ok <- is.finite(x) & is.finite(y)
  x <- x[ok]; y <- y[ok]
  ord <- order(x); x <- x[ord]; y <- y[ord]
  n <- length(x)
  if (n < 2L * min_n + 2L) return(NULL)

  ux <- sort(unique(x))
  if (is.null(grid)) {
    lo <- x[min_n]
    hi <- x[n - min_n + 1L]
    grid <- ux[ux > lo & ux < hi]
  }
  if (!length(grid)) return(NULL)

  null_fit <- try(
    if (outcome_type == "binary") stats::glm(y ~ x, family = stats::binomial())
    else if (outcome_type == "count") stats::glm(y ~ x, family = stats::poisson())
    else stats::lm(y ~ x),
    silent = TRUE
  )
  null_bic <- if (inherits(null_fit, "try-error")) NA_real_ else stats::BIC(null_fit)

  fits <- lapply(grid, function(bp) .tablearn_fit_hinge(x, y, bp, outcome_type))
  bics <- vapply(fits, function(f) if (is.null(f)) Inf else stats::BIC(f), numeric(1))
  if (!any(is.finite(bics))) return(NULL)
  j <- which.min(bics)
  fit <- fits[[j]]; bp <- as.numeric(grid[j])
  cf <- stats::coef(fit)
  slope1 <- unname(cf["x"])
  slope_change <- unname(cf["hinge"])
  slope2 <- slope1 + slope_change
  sm <- summary(fit)$coefficients
  p_change <- if ("hinge" %in% rownames(sm)) sm["hinge", ncol(sm)] else NA_real_

  out <- list(
    breakpoint = bp,
    bic = bics[j],
    bic_no_break = null_bic,
    delta_bic = null_bic - bics[j],
    detected = is.finite(null_bic) && (bics[j] + 2 < null_bic),
    slope_before = slope1,
    slope_after = slope2,
    slope_change = slope_change,
    p_change = p_change,
    model = fit,
    ci = c(NA_real_, NA_real_),
    boot = numeric()
  )

  if (is.numeric(boot) && boot > 0L) {
    bps <- numeric()
    for (b in seq_len(as.integer(boot))) {
      ii <- sample.int(n, n, replace = TRUE)
      xb <- x[ii]; yb <- y[ii]
      oo <- order(xb); xb <- xb[oo]; yb <- yb[oo]
      tmp <- .tablearn_breakpoint(xb, yb, outcome_type, min_n, grid = grid, boot = 0L, ci_level = ci_level)
      if (!is.null(tmp) && is.finite(tmp$breakpoint)) bps <- c(bps, tmp$breakpoint)
    }
    if (length(bps) >= 20L) {
      a <- (1 - ci_level) / 2
      out$ci <- as.numeric(stats::quantile(bps, c(a, 1 - a), na.rm = TRUE, names = FALSE))
      out$boot <- bps
    }
  }
  out
}

.tablearn_cusum <- function(y, outcome_type, target, alt = NULL, better = "lower",
                            method = "deviation", reset = FALSE, limit = NULL) {
  if (is.null(target) || !is.numeric(target) || length(target) != 1L || !is.finite(target)) {
    .tablearn_stop("`target` must be supplied for CUSUM.")
  }
  y <- as.numeric(y)

  if (method == "llr") {
    if (outcome_type != "binary") .tablearn_stop("Likelihood-ratio CUSUM is currently implemented for binary outcomes only.")
    if (is.null(alt) || !is.numeric(alt) || length(alt) != 1L || alt <= 0 || alt >= 1) {
      .tablearn_stop("For binary `cusum_method = 'llr'`, supply `cusum_alt` in (0, 1).")
    }
    if (target <= 0 || target >= 1) .tablearn_stop("Binary CUSUM target must lie in (0, 1).")
    score <- ifelse(y == 1,
      log(alt / target),
      log((1 - alt) / (1 - target))
    )
  } else {
    score <- if (better == "higher") target - y else y - target
  }

  cs <- numeric(length(score))
  prev <- 0
  for (i in seq_along(score)) {
    if (is.na(score[i])) {
      cs[i] <- prev
      next
    }
    val <- prev + score[i]
    if (reset) val <- max(0, val)
    prev <- val
    cs[i] <- val
  }
  signal <- if (!is.null(limit) && is.finite(limit)) which(cs >= limit)[1L] else NA_integer_
  data.frame(.case = seq_along(y), .value = cs, .score = score, .signal = seq_along(y) == signal)
}

.tablearn_lccusum <- function(y_event, p_acceptable, p_unacceptable, alpha = 0.05, beta = 0.10) {
  if (is.null(p_acceptable) || is.null(p_unacceptable)) {
    .tablearn_stop("`p_acceptable` and `p_unacceptable` are required for LC-CUSUM.")
  }
  if (p_acceptable <= 0 || p_acceptable >= 1 || p_unacceptable <= 0 || p_unacceptable >= 1) {
    .tablearn_stop("LC-CUSUM failure rates must lie in (0, 1).")
  }
  if (p_acceptable >= p_unacceptable) {
    .tablearn_stop("For LC-CUSUM based on failure rates, `p_acceptable` must be smaller than `p_unacceptable`.")
  }
  P <- log(p_acceptable / p_unacceptable)
  Q <- log((1 - p_unacceptable) / (1 - p_acceptable))
  s <- Q / (P + Q)
  a <- log((1 - beta) / alpha)
  b <- log((1 - alpha) / beta)
  H <- a / (P + Q) # negative lower decision limit in this parameterization
  score <- as.numeric(y_event) - s
  cs <- cumsum(ifelse(is.na(score), 0, score))
  prof <- which(cs <= H)[1L]
  data.frame(
    .case = seq_along(y_event), .value = cs, .score = score,
    .limit = H, .proficient = seq_along(y_event) == prof
  ) -> dat
  attr(dat, "parameters") <- list(P = P, Q = Q, s = s, a = a, b = b, H = H, proficiency_case = prof)
  dat
}

.tablearn_racusum <- function(y_event, p_expected, odds_ratio = 2, limit = NULL) {
  if (length(y_event) != length(p_expected)) .tablearn_stop("`expected` risk must have one value per analyzed case.")
  if (any(p_expected <= 0 | p_expected >= 1, na.rm = TRUE)) {
    .tablearn_stop("Expected risks for RA-CUSUM must lie strictly between 0 and 1.")
  }
  if (!is.numeric(odds_ratio) || length(odds_ratio) != 1L || odds_ratio <= 1) {
    .tablearn_stop("`racusum_or` must be > 1 when monitoring deterioration.")
  }
  y <- as.numeric(y_event)
  p <- as.numeric(p_expected)
  R0 <- 1
  RA <- odds_ratio
  w_fail <- log(((1 - p + R0 * p) * RA) / ((1 - p + RA * p) * R0))
  w_success <- log((1 - p + R0 * p) / (1 - p + RA * p))
  score <- ifelse(y == 1, w_fail, w_success)
  cs <- numeric(length(score)); prev <- 0
  for (i in seq_along(score)) {
    if (is.na(score[i])) { cs[i] <- prev; next }
    prev <- max(0, prev + score[i])
    cs[i] <- prev
  }
  sig <- if (!is.null(limit) && is.finite(limit)) which(cs >= limit)[1L] else NA_integer_
  data.frame(.case = seq_along(y), .value = cs, .score = score, .expected = p, .signal = seq_along(y) == sig)
}

.tablearn_phase_table <- function(raw, breaks, outcome_type, phase_names = NULL, ci_level = 0.95) {
  if (is.null(breaks) || !length(breaks) || !nrow(raw)) return(data.frame())
  breaks <- sort(unique(as.numeric(breaks[is.finite(breaks)])))
  cuts <- c(-Inf, breaks, Inf)
  k <- length(cuts) - 1L
  if (is.null(phase_names)) phase_names <- c("Learning", "Consolidation", "Proficiency")
  if (k == 2L && length(phase_names) >= 3L) {
    phase_names <- c(phase_names[1L], phase_names[3L])
  } else if (length(phase_names) < k) {
    phase_names <- c(phase_names, paste0("Phase ", seq_len(k - length(phase_names))))
  }
  phase_names <- phase_names[seq_len(k)]
  ph <- cut(raw$.case, breaks = cuts, labels = phase_names, right = TRUE)
  spl <- split(seq_len(nrow(raw)), ph, drop = TRUE)
  do.call(rbind, lapply(names(spl), function(nm) {
    ii <- spl[[nm]]; y <- raw$.value[ii]; n <- sum(!is.na(y))
    if (outcome_type == "binary") {
      ev <- sum(y, na.rm = TRUE); est <- mean(y, na.rm = TRUE); ci <- .tablearn_wilson(ev, n, ci_level)
      data.frame(phase = nm, start = min(raw$.case[ii]), end = max(raw$.case[ii]), n = n,
                 events = ev, estimate = est, lower = ci[1], upper = ci[2])
    } else {
      est <- mean(y, na.rm = TRUE); sdv <- stats::sd(y, na.rm = TRUE)
      q <- stats::qt(1 - (1 - ci_level) / 2, df = max(1, n - 1)); se <- sdv / sqrt(n)
      data.frame(phase = nm, start = min(raw$.case[ii]), end = max(raw$.case[ii]), n = n,
                 events = NA_real_, estimate = est, lower = est - q * se, upper = est + q * se)
    }
  }))
}


.tablearn_first_sustained <- function(ok, cases, hold = 3L) {
  hold <- as.integer(hold)
  if (length(hold) != 1L || is.na(hold) || hold < 1L) {
    .tablearn_stop("`target_hold` must be one integer >= 1.")
  }
  ok <- as.logical(ok)
  ok[is.na(ok)] <- FALSE
  cases <- as.numeric(cases)
  if (!length(ok) || length(ok) != length(cases)) {
    return(list(first = NA_real_, confirmed = NA_real_, current = FALSE, current_run = 0L))
  }
  rr <- rle(ok)
  ends <- cumsum(rr$lengths)
  starts <- ends - rr$lengths + 1L
  jj <- which(rr$values & rr$lengths >= hold)[1L]
  first <- confirmed <- NA_real_
  if (!is.na(jj)) {
    first <- cases[starts[jj]]
    confirmed <- cases[starts[jj] + hold - 1L]
  }
  current_run <- if (length(rr$values) && isTRUE(tail(rr$values, 1L))) tail(rr$lengths, 1L) else 0L
  list(
    first = first,
    confirmed = confirmed,
    current = current_run >= hold,
    current_run = as.integer(current_run)
  )
}

.tablearn_variability <- function(y, metric = "sd") {
  y <- as.numeric(y)
  y <- y[is.finite(y)]
  if (length(y) < 2L) return(NA_real_)
  if (metric == "sd") return(stats::sd(y))
  if (metric == "iqr") return(stats::IQR(y, na.rm = TRUE, type = 7))
  if (metric == "cv") {
    mm <- mean(y)
    if (!is.finite(mm) || abs(mm) < sqrt(.Machine$double.eps)) return(NA_real_)
    return(stats::sd(y) / abs(mm))
  }
  NA_real_
}

.tablearn_proficiency_engine <- function(
    raw, win, bp, lccusum, racusum,
    outcome_type, better, target,
    proficiency = TRUE,
    proficiency_method = "auto",
    proficiency_case = NULL,
    proficiency_case_rule = "confirmed",
    target_hold = 3L,
    target_tolerance = 0,
    plateau = TRUE,
    plateau_ratio = 0.25,
    plateau_slope = NULL,
    stability = TRUE,
    stability_metric = "auto",
    stability_window = NULL,
    stability_ratio = 0.75,
    window_n = 5L,
    operator = "Overall") {

  if (!isTRUE(proficiency)) {
    sm <- data.frame(
      operator = operator, status = "Not assessed", achieved = FALSE,
      proficiency_case = NA_real_, method = proficiency_method,
      evidence = "Not assessed", breakpoint_case = NA_real_,
      breakpoint_detected = FALSE, plateau = NA, slope_before = NA_real_,
      slope_after = NA_real_, target = if (is.null(target)) NA_real_ else target,
      target_first_case = NA_real_, target_confirmed_case = NA_real_,
      target_achieved = NA, target_current = NA,
      lccusum_case = NA_real_, lccusum_achieved = NA,
      stability = NA, stability_metric = NA_character_,
      variability_early = NA_real_, variability_late = NA_real_,
      variability_ratio = NA_real_, current_case = if (nrow(raw)) max(raw$.case) else NA_real_,
      current_phase = "Not assessed", deterioration_signal = NA,
      deterioration_case = NA_real_, cases_since_proficiency = NA_real_,
      stringsAsFactors = FALSE
    )
    return(list(summary = sm, evidence = data.frame(), phase_breaks = numeric(), phase_names = character(), current_status = data.frame()))
  }

  proficiency_method <- match.arg(
    proficiency_method,
    c("auto", "combined", "target", "lccusum", "breakpoint", "manual")
  )
  proficiency_case_rule <- match.arg(proficiency_case_rule, c("confirmed", "first_stable"))
  stability_metric <- match.arg(stability_metric, c("auto", "sd", "iqr", "cv", "none"))

  if (!is.numeric(target_tolerance) || length(target_tolerance) != 1L || is.na(target_tolerance) || target_tolerance < 0) {
    .tablearn_stop("`target_tolerance` must be one non-negative number.")
  }
  if (!is.numeric(plateau_ratio) || length(plateau_ratio) != 1L || is.na(plateau_ratio) || plateau_ratio < 0 || plateau_ratio > 1) {
    .tablearn_stop("`plateau_ratio` must lie in [0, 1].")
  }
  if (!is.null(plateau_slope) && (!is.numeric(plateau_slope) || length(plateau_slope) != 1L || !is.finite(plateau_slope) || plateau_slope < 0)) {
    .tablearn_stop("`plateau_slope` must be NULL or one non-negative finite number.")
  }
  if (!is.numeric(stability_ratio) || length(stability_ratio) != 1L || is.na(stability_ratio) || stability_ratio <= 0) {
    .tablearn_stop("`stability_ratio` must be one positive number.")
  }

  current_case <- if (nrow(raw)) max(raw$.case, na.rm = TRUE) else NA_real_

  # Breakpoint / plateau evidence. A change point alone is not called proficiency.
  break_detected <- !is.null(bp) && isTRUE(bp$detected) && is.finite(bp$breakpoint)
  break_case <- if (break_detected) as.numeric(bp$breakpoint) else NA_real_
  slope_before <- if (!is.null(bp) && is.finite(bp$slope_before)) as.numeric(bp$slope_before) else NA_real_
  slope_after <- if (!is.null(bp) && is.finite(bp$slope_after)) as.numeric(bp$slope_after) else NA_real_
  plateau_met <- FALSE
  plateau_threshold <- NA_real_
  if (isTRUE(plateau) && break_detected && is.finite(slope_before) && is.finite(slope_after)) {
    learning_direction <- if (better == "higher") slope_before > 0 else slope_before < 0
    plateau_threshold <- if (is.null(plateau_slope)) abs(slope_before) * plateau_ratio else plateau_slope
    plateau_met <- isTRUE(learning_direction) && is.finite(plateau_threshold) && abs(slope_after) <= plateau_threshold
  }
  plateau_case <- if (plateau_met) break_case else NA_real_

  # Sustained target on the rolling/block/current-performance series.
  target_available <- !is.null(target) && is.numeric(target) && length(target) == 1L && is.finite(target)
  target_run <- list(first = NA_real_, confirmed = NA_real_, current = FALSE, current_run = 0L)
  if (target_available && nrow(win)) {
    target_ok <- if (better == "higher") {
      win$.value >= (target - target_tolerance)
    } else {
      win$.value <= (target + target_tolerance)
    }
    target_run <- .tablearn_first_sustained(target_ok, win$.end, target_hold)
  }
  target_case <- if (proficiency_case_rule == "first_stable") target_run$first else target_run$confirmed
  target_achieved <- target_available && is.finite(target_case)

  # LC-CUSUM is a direct sequential competency criterion for binary outcomes.
  lc_case <- NA_real_
  if (is.data.frame(lccusum) && nrow(lccusum)) {
    jj <- which(lccusum$.proficient %in% TRUE)[1L]
    if (!is.na(jj)) lc_case <- as.numeric(lccusum$.case[jj])
  }
  lc_available <- is.data.frame(lccusum) && nrow(lccusum) > 0L
  lc_achieved <- lc_available && is.finite(lc_case)

  # Stability compares equally sized early and late raw-case segments. For binary
  # outcomes it is off by default because SD is largely determined by the event rate.
  metric_resolved <- stability_metric
  if (metric_resolved == "auto") metric_resolved <- if (outcome_type == "binary") "none" else "sd"
  stability_met <- NA
  var_early <- var_late <- var_ratio <- NA_real_
  if (isTRUE(stability) && metric_resolved != "none" && nrow(raw) >= 6L) {
    sw <- if (is.null(stability_window)) max(5L, as.integer(window_n)) else as.integer(stability_window)
    if (length(sw) != 1L || is.na(sw) || sw < 3L) .tablearn_stop("`stability_window` must be NULL or one integer >= 3.")
    sw <- min(sw, floor(nrow(raw) / 2L))
    if (sw >= 3L) {
      var_early <- .tablearn_variability(head(raw$.value, sw), metric_resolved)
      var_late <- .tablearn_variability(tail(raw$.value, sw), metric_resolved)
      if (is.finite(var_early) && is.finite(var_late)) {
        if (var_early <= sqrt(.Machine$double.eps)) {
          var_ratio <- if (var_late <= sqrt(.Machine$double.eps)) 1 else Inf
        } else {
          var_ratio <- var_late / var_early
        }
        stability_met <- is.finite(var_ratio) && var_ratio <= stability_ratio
      }
    }
  }

  # A deterioration signal does not erase prior proficiency, but it changes the
  # current status and should be visible to the user.
  deterioration_case <- NA_real_
  if (is.data.frame(racusum) && nrow(racusum)) {
    jj <- which(racusum$.signal %in% TRUE)[1L]
    if (!is.na(jj)) deterioration_case <- as.numeric(racusum$.case[jj])
  }
  deterioration_signal <- is.finite(deterioration_case)

  # Resolve the requested proficiency rule.
  method_used <- if (proficiency_method == "auto") "combined" else proficiency_method
  candidate <- NA_real_
  status <- "Insufficient evidence"

  if (method_used == "manual") {
    if (is.null(proficiency_case) || !is.numeric(proficiency_case) || length(proficiency_case) != 1L || !is.finite(proficiency_case)) {
      .tablearn_stop("`proficiency_case` must be supplied as one finite case number when `proficiency_method = 'manual'`.")
    }
    candidate <- as.numeric(proficiency_case)
    status <- "User-defined"
  } else if (method_used == "target") {
    if (!target_available) .tablearn_stop("`target` is required when `proficiency_method = 'target'`.")
    if (target_achieved) { candidate <- target_case; status <- "Achieved" } else status <- "Not achieved"
  } else if (method_used == "lccusum") {
    if (!lc_available) .tablearn_stop("Enable `lccusum = TRUE` when `proficiency_method = 'lccusum'`.")
    if (lc_achieved) { candidate <- lc_case; status <- "Achieved" } else status <- "Not achieved"
  } else if (method_used == "breakpoint") {
    cc <- if (isTRUE(plateau)) plateau_case else break_case
    if (is.finite(cc)) { candidate <- cc; status <- "Achieved" }
    else if (break_detected) status <- "Not achieved"
  } else {
    # Combined/auto: explicitly requested clinical/sequential criteria are treated
    # as requirements. Plateau and stability strengthen evidence but are not allowed
    # to override failure to meet a supplied target or LC-CUSUM criterion.
    required_missing <- FALSE
    cases <- numeric()
    if (target_available) {
      if (target_achieved) cases <- c(cases, target_case) else required_missing <- TRUE
    }
    if (lc_available) {
      if (lc_achieved) cases <- c(cases, lc_case) else required_missing <- TRUE
    }
    if (required_missing) {
      status <- "Not achieved"
    } else {
      if (plateau_met) cases <- c(cases, plateau_case)
      if (length(cases)) {
        candidate <- max(cases, na.rm = TRUE)
        status <- "Achieved"
      } else if (plateau_met) {
        candidate <- plateau_case
        status <- "Achieved"
      }
    }
  }

  achieved <- is.finite(candidate) && status %in% c("Achieved", "User-defined")

  evidence_met <- c(
    sustained_target = isTRUE(target_achieved),
    lccusum = isTRUE(lc_achieved),
    plateau = isTRUE(plateau_met),
    stability = isTRUE(stability_met)
  )
  n_evidence <- sum(evidence_met)
  evidence_level <- if (status == "User-defined") {
    "User-defined"
  } else if (!achieved) {
    if (status == "Not achieved") "Not achieved" else "Insufficient"
  } else if (n_evidence >= 3L) {
    "Strong"
  } else if (n_evidence == 2L) {
    "Moderate"
  } else {
    "Limited"
  }

  current_phase <- "Learning"
  if (achieved && is.finite(current_case) && current_case >= candidate) {
    current_phase <- "Proficiency"
  } else if (break_detected && is.finite(current_case) && current_case > break_case) {
    current_phase <- "Consolidation"
  }
  if (deterioration_signal && is.finite(current_case) && deterioration_case <= current_case) {
    current_phase <- "Deterioration signal"
  }

  # Automatic phase boundaries are descriptive labels, not additional hypothesis
  # tests. The proficiency boundary starts at the estimated proficiency case.
  auto_breaks <- numeric()
  auto_names <- character()
  if (achieved) {
    prof_boundary <- candidate - 1
    if (break_detected && is.finite(break_case) && break_case < prof_boundary) {
      auto_breaks <- c(break_case, prof_boundary)
      auto_names <- c("Learning", "Consolidation", "Proficiency")
    } else if (is.finite(prof_boundary) && prof_boundary >= min(raw$.case, na.rm = TRUE)) {
      auto_breaks <- prof_boundary
      auto_names <- c("Learning", "Proficiency")
    }
  } else if (break_detected) {
    auto_breaks <- break_case
    auto_names <- c("Learning", "Consolidation")
  }

  evidence <- data.frame(
    operator = operator,
    criterion = c("Change point", "Plateau", "Sustained target", "LC-CUSUM", "Stability"),
    available = c(!is.null(bp), isTRUE(plateau) && break_detected, target_available, lc_available,
                  isTRUE(stability) && metric_resolved != "none"),
    met = c(isTRUE(break_detected), isTRUE(plateau_met), isTRUE(target_achieved), isTRUE(lc_achieved), isTRUE(stability_met)),
    case = c(break_case, plateau_case, target_case, lc_case, NA_real_),
    detail = c(
      if (break_detected) paste0("Change point near case ", signif(break_case, 4)) else "No clear change point",
      if (isTRUE(plateau)) paste0("Post-change |slope| threshold = ", signif(plateau_threshold, 4)) else "Not assessed",
      if (target_available) paste0("Target ", signif(target, 4), "; hold = ", as.integer(target_hold), " windows") else "No target supplied",
      if (lc_available) "LC-CUSUM enabled" else "LC-CUSUM not enabled",
      if (metric_resolved != "none") paste0(metric_resolved, " late/early ratio = ", signif(var_ratio, 4)) else "Not assessed"
    ),
    stringsAsFactors = FALSE
  )

  summary <- data.frame(
    operator = operator,
    status = status,
    achieved = achieved,
    proficiency_case = if (achieved) candidate else NA_real_,
    method = method_used,
    evidence = evidence_level,
    breakpoint_case = break_case,
    breakpoint_detected = break_detected,
    plateau = plateau_met,
    slope_before = slope_before,
    slope_after = slope_after,
    target = if (target_available) target else NA_real_,
    target_first_case = target_run$first,
    target_confirmed_case = target_run$confirmed,
    target_achieved = if (target_available) target_achieved else NA,
    target_current = if (target_available) target_run$current else NA,
    lccusum_case = lc_case,
    lccusum_achieved = if (lc_available) lc_achieved else NA,
    stability = stability_met,
    stability_metric = metric_resolved,
    variability_early = var_early,
    variability_late = var_late,
    variability_ratio = var_ratio,
    current_case = current_case,
    current_phase = current_phase,
    deterioration_signal = deterioration_signal,
    deterioration_case = deterioration_case,
    cases_since_proficiency = if (achieved && is.finite(current_case)) max(0, current_case - candidate) else NA_real_,
    stringsAsFactors = FALSE
  )

  current_status <- summary[, c(
    "operator", "current_case", "current_phase", "proficiency_case",
    "cases_since_proficiency", "target_current", "deterioration_signal", "deterioration_case"
  ), drop = FALSE]

  list(
    summary = summary,
    evidence = evidence,
    phase_breaks = auto_breaks,
    phase_names = auto_names,
    current_status = current_status
  )
}

.tablearn_recursive_modify <- function(base, extra) {
  if (is.null(extra)) return(base)
  for (nm in names(extra)) {
    if (is.list(extra[[nm]]) && is.list(base[[nm]])) base[[nm]] <- .tablearn_recursive_modify(base[[nm]], extra[[nm]])
    else base[[nm]] <- extra[[nm]]
  }
  base
}

.tablearn_plot_defaults <- function() {
  list(
    raw = list(show = TRUE, color = "#7A7A7A", fill = NA, shape = 16, size = 1.5, alpha = 0.28, stroke = 0.3),
    window = list(show = TRUE, geom = "point_line", color = "#1F5A94", fill = "white", shape = 21,
                  size = 2.6, alpha = 1, line_color = "#1F5A94", line_width = 1.05, line_type = 1),
    window_interval = list(show = TRUE, geom = "ribbon", color = "#1F5A94", fill = "#1F5A94",
                           alpha = 0.13, line_width = 0.45, width = 0.25),
    smooth = list(show = TRUE, color = "#B23A48", fill = "#B23A48", line_width = 1.2, line_type = 1,
                  alpha = 1, ci_show = TRUE, ci_alpha = 0.10),
    ewma = list(show = TRUE, color = "#7A3E9D", line_width = 1.1, line_type = 2, alpha = 1),
    breakpoint = list(show = TRUE, color = "#222222", line_width = 0.8, line_type = 2, alpha = 0.9,
                      label = TRUE, label_text = NULL, label_size = 3.4, label_angle = 90,
                      label_hjust = -0.08, label_vjust = 1.15),
    proficiency = list(show = TRUE, color = "#00796B", line_width = 1.0, line_type = 1, alpha = 0.95,
                       label = TRUE, label_text = NULL, label_size = 3.4, label_angle = 90,
                       label_hjust = -0.08, label_vjust = 1.15),
    phase = list(show = TRUE, fills = c("#F5B7B1", "#F9E79F", "#ABEBC6"), alpha = 0.11,
                 border_color = NA, border_width = 0, label = TRUE, label_size = 3.4,
                 label_position = "top"),
    target = list(show = TRUE, color = "#2E7D32", line_width = 0.9, line_type = 3, alpha = 0.9,
                  label = TRUE, label_text = NULL, label_size = 3.3),
    current = list(show = TRUE, color = "#000000", fill = "#FFD166", shape = 21, size = 4,
                   alpha = 1, stroke = 0.8, label = TRUE, label_text = NULL, label_size = 3.4,
                   hjust = -0.08, vjust = -0.7),
    axes = list(x_limits = NULL, y_limits = NULL, x_breaks = NULL, y_breaks = NULL,
                x_expand = NULL, y_expand = NULL, y_percent = "auto", percent_accuracy = 1,
                x_reverse = FALSE, y_reverse = FALSE, x_trans = NULL, y_trans = NULL,
                x_date_format = NULL, x_date_breaks = NULL, clip = "on"),
    legend = list(show = TRUE, position = "bottom", title = NULL, direction = "horizontal"),
    facet = list(ncol = NULL, nrow = NULL, scales = "fixed"),
    operator = list(colors = NULL, shapes = NULL, line_types = NULL),
    cusum = list(
      colors = c("CUSUM" = "#1F5A94", "LC-CUSUM" = "#B23A48", "RA-CUSUM" = "#7A3E9D"),
      line_types = c("CUSUM" = 1, "LC-CUSUM" = 1, "RA-CUSUM" = 1),
      line_width = 1.05, alpha = 1, points = FALSE, point_size = 1.5, point_shape = 16,
      zero_show = TRUE, zero_color = "#777777", zero_width = 0.45, zero_type = 1,
      decision_show = TRUE, decision_color = "#222222", decision_width = 0.75, decision_type = 2,
      signal_show = TRUE, signal_color = "#000000", signal_fill = "#FFD166",
      signal_shape = 21, signal_size = 3.8, signal_stroke = 0.8,
      title = NULL, subtitle = NULL, caption = NULL, xlab = NULL, ylab = "CUSUM",
      x_limits = NULL, y_limits = NULL, x_breaks = NULL, y_breaks = NULL
    ),
    theme = list(name = "publication", base_size = 11, base_family = "", grid_major = TRUE,
                 grid_minor = FALSE, panel_border = TRUE, axis_line = FALSE,
                 plot_title_face = "bold", legend_key_size = 0.9),
    text = list(title_size = NULL, subtitle_size = NULL, caption_size = NULL,
                axis_title_size = NULL, axis_text_size = NULL, legend_text_size = NULL),
    margins = list(top = 8, right = 12, bottom = 8, left = 8),
    panel = list(background = NULL, border_color = "#333333", border_width = 0.5),
    annotation = list(show_n = FALSE, show_window_label = FALSE)
  )
}

#' Learning-curve analysis for sequential clinical or procedural performance
#'
#' `tablearn()` analyzes performance across consecutive cases, procedures, or
#' observations. It keeps individual case-level data for statistical inference
#' while allowing visually stable rolling or block summaries for learning-curve
#' display. A right-aligned rolling window such as cases 1-5, 2-6, 3-7, ... is
#' particularly useful when the question is: "What is my current performance
#' based on the most recent N cases?"
#'
#' @param outcome Outcome variable. May be a bare variable name, character name,
#'   or `vars(...)` containing multiple outcomes.
#' @param data Data frame. If `NULL`, `tablearn()` attempts to use the R4VN active
#'   data frame.
#' @param order Optional ordering variable such as case number or procedure date.
#'   If omitted, current row order is used. Data are sorted by this variable within
#'   each operator.
#' @param operator Optional operator/surgeon/trainee variable. Curves and analyses
#'   are then calculated separately within operator.
#' @param type Outcome type: `"auto"`, `"continuous"`, `"binary"`, or `"count"`.
#'   In auto mode, logical/factor/character variables with exactly two non-missing
#'   levels and numeric 0/1 variables are treated as binary. A numeric variable
#'   with two other observed values (for example, 50 and 100) is treated as
#'   continuous unless `event` is supplied or `type = "binary"` is requested.
#'   This conservative rule prevents a two-valued continuous measurement from
#'   being silently converted to a proportion.
#' @param event Event level for a binary outcome. If omitted, the second factor
#'   level, `TRUE`, or the larger numeric level is used.
#' @param better Direction of better performance: `"lower"`, `"higher"`, or
#'   `"auto"`. Auto uses `"lower"` for numeric outcomes and for an event coded as
#'   failure/complication it should normally be set explicitly by the analyst.
#' @param window Number of cases in a rolling or block window. Default 5.
#' @param window_type `"rolling"` (1-5, 2-6, 3-7, ...), `"block"` (1-5,
#'   6-10, ...), `"cumulative"` (1, 1-2, 1-3, ...), or `"none"`.
#' @param window_align Alignment for rolling windows: `"right"` (recommended and
#'   default for current-performance monitoring), `"center"`, or `"left"`.
#' @param window_step Number of cases to advance each window. Default is 1 for
#'   rolling/cumulative and `window` for block summaries.
#' @param window_complete If `TRUE`, only complete windows are retained. If
#'   `FALSE`, partial windows at the beginning/end are allowed.
#' @param window_fun `"auto"`, `"mean"`, `"median"`, `"trimmed_mean"`,
#'   `"proportion"`, or `"rate"`. Auto uses a proportion for binary outcomes and
#'   mean otherwise.
#' @param window_trim Trim proportion used by `window_fun = "trimmed_mean"`.
#' @param window_weight `"equal"`, `"linear"`, or `"exponential"`. Equal weighting
#'   gives the ordinary rolling mean/proportion; alternatives give more weight to
#'   recent cases within each window.
#' @param window_decay Decay in `(0,1]` for exponentially weighted windows.
#' @param window_ci Show point-wise interval for each window.
#' @param ci_level Confidence level, default 0.95.
#' @param window_ci_method `"auto"`, `"t"`, `"normal"`, `"wilson"`, or `"exact"`.
#'   Auto uses Wilson for binary proportions, t intervals for small continuous
#'   windows, and normal intervals for larger continuous windows.
#' @param interval Window interval displayed/calculated: `"ci"`, `"sd"`, `"iqr"`,
#'   `"range"`, or `"none"`.
#' @param ewma Logical; calculate an exponentially weighted moving average.
#' @param ewma_lambda EWMA smoothing parameter in `(0,1]`; larger values react more
#'   strongly to the latest case.
#' @param ewma_init EWMA starting value: `"first"`, `"mean"`, or `"target"`.
#' @param smooth Logical; add a fitted smooth trend.
#' @param smooth_method `"loess"`, `"spline"`, `"lm"`, `"gam"`, or `"none"`.
#'   GAM requires package `mgcv`; otherwise loess is used as fallback.
#' @param smooth_on Data used only for the descriptive smooth: `"window"`, `"raw"`,
#'   or `"ewma"`.
#' @param smooth_span LOESS span.
#' @param smooth_df Optional degrees of freedom for smoothing spline.
#' @param smooth_ci Show a point-wise confidence band where the smoothing method
#'   supports one.
#' @param breakpoint Logical; estimate a one-change-point piecewise regression.
#' @param break_on `"raw"` (recommended default for inference) or `"window"`.
#'   Overlapping rolling windows are correlated, so fitting the inferential model
#'   to raw cases avoids treating overlapping windows as independent observations.
#' @param break_min_n Minimum observations required on each side of a candidate
#'   breakpoint.
#' @param break_grid Optional numeric vector of candidate case positions.
#' @param break_boot Number of bootstrap replications for an exploratory percentile
#'   CI for the breakpoint. `0` disables bootstrap.
#' @param phase Logical; create a phase summary table using the estimated or user
#'   supplied phase breaks.
#' @param phase_breaks Optional numeric vector of manual phase boundaries. This is
#'   useful for 3+ named phases; automatic estimation currently provides one main
#'   change point.
#' @param phase_names Names of phases. Defaults include Learning, Consolidation,
#'   and Proficiency. Manual names are used when `phase_breaks` is supplied.
#' @param phase_method Phase classification rule: `"auto"`/`"combined"` uses the
#'   detected learning change point together with the estimated proficiency case;
#'   `"manual"` uses only `phase_breaks`; `"breakpoint"` uses the change point;
#'   and `"proficiency"` divides the series at the proficiency case. Automatic
#'   classification never labels a phase Proficiency unless proficiency is actually
#'   estimated.
#' @param proficiency Logical; assess whether proficiency has been reached.
#' @param proficiency_method `"auto"`, `"combined"`, `"target"`, `"lccusum"`,
#'   `"breakpoint"`, or `"manual"`. The recommended `"auto"` resolves to a
#'   combined rule. A supplied target and an enabled LC-CUSUM are treated as
#'   required criteria; plateau/stability strengthen the evidence.
#' @param proficiency_case Manual case number used only with
#'   `proficiency_method = "manual"`.
#' @param proficiency_case_rule For sustained-target proficiency, report the first
#'   qualifying window (`"first_stable"`) or the window in which the required run
#'   is confirmed (`"confirmed"`, default).
#' @param target_hold Number of consecutive window estimates that must satisfy
#'   `target` before target-based proficiency is confirmed. Default 3.
#' @param target_tolerance Non-negative tolerance around `target`. For
#'   `better = "lower"`, values <= target + tolerance qualify; for
#'   `better = "higher"`, values >= target - tolerance qualify.
#' @param plateau Logical; assess whether the post-breakpoint slope is sufficiently
#'   small to be interpreted as an operational plateau. A breakpoint alone is not
#'   automatically called proficiency.
#' @param plateau_ratio Relative plateau threshold when `plateau_slope` is `NULL`.
#'   The post-breakpoint absolute slope must be <= `plateau_ratio` times the initial
#'   absolute slope. Default 0.25.
#' @param plateau_slope Optional absolute slope threshold overriding
#'   `plateau_ratio`. This is outcome-scale specific and can be useful when a
#'   clinically meaningful slope threshold is known.
#' @param stability Logical; assess whether performance variability has fallen.
#' @param stability_metric `"auto"`, `"sd"`, `"iqr"`, `"cv"`, or `"none"`.
#'   Auto uses SD for continuous/count outcomes and does not use variability as
#'   proficiency evidence for binary outcomes.
#' @param stability_window Number of early and late raw cases used to compare
#'   variability. Default is at least the selected learning-curve window size.
#' @param stability_ratio Late/early variability ratio required for stability.
#'   Default 0.75 means late variability must be at most 75% of early variability.
#' @param cusum Logical; calculate a conventional CUSUM from individual cases.
#' @param target Clinical/quality target for CUSUM, target/reference line, and
#'   optional EWMA initialization.
#' @param cusum_method `"deviation"` or binary likelihood-ratio `"llr"`.
#' @param cusum_alt Alternative binary failure probability for LLR-CUSUM.
#' @param cusum_reset If `TRUE`, use a one-sided tabular CUSUM reset at zero.
#' @param cusum_limit Optional decision limit. If omitted, CUSUM is descriptive.
#' @param lccusum Logical; calculate a binary LC-CUSUM designed to signal evidence
#'   that an acceptable failure rate has been reached.
#' @param p_acceptable Acceptable failure probability for LC-CUSUM.
#' @param p_unacceptable Unacceptable failure probability for LC-CUSUM; must be
#'   larger than `p_acceptable`.
#' @param alpha,beta Type-I and Type-II error probabilities used in the LC-CUSUM
#'   decision boundary formula.
#' @param racusum Logical; calculate a binary risk-adjusted CUSUM using likelihood
#'   scores and individual expected risks.
#' @param expected Optional expected-risk variable or numeric vector for RA-CUSUM.
#'   Supplying externally validated or pre-operative expected risks is preferable.
#' @param adjust Optional covariates used to fit a logistic expected-risk model if
#'   `expected` is not supplied. May be character names or `vars(...)`.
#' @param racusum_or Odds ratio representing deterioration to be detected. Must be
#'   greater than 1.
#' @param racusum_limit Optional RA-CUSUM decision limit.
#' @param sensitivity Logical; calculate window-size sensitivity summaries.
#' @param window_sensitivity Numeric vector of window sizes. If omitted and
#'   `sensitivity = TRUE`, sensible values are chosen from the sample size.
#' @param plot Logical; create figures. Standard figures and Viewer figures use
#'   base R and therefore require no add-on package. When `ggplot2` is installed,
#'   an advanced ggplot object is also retained in `$plots`.
#' @param plot_type `"auto"`, `"learning"`, `"cusum"`, or `"both"`. Auto shows
#'   the learning curve and also the CUSUM panel whenever conventional CUSUM,
#'   LC-CUSUM, or RA-CUSUM has been requested.
#' @param x_axis `"case"` (default) or `"order"`. Case number is usually preferable
#'   for learning curves; date/time can be shown using `"order"`.
#' @param operator_display `"facet"` or `"overlay"` when `operator` is supplied.
#' @param show_raw,show_window,show_smooth,show_ewma,show_interval,show_break,show_proficiency,show_phase,show_target,show_current High-level plot layer switches.
#' @param title,subtitle,caption,xlab,ylab Plot labels. Defaults are generated from
#'   the outcome and selected language.
#' @param theme Plot theme: `"publication"`, `"minimal"`, `"classic"`, `"bw"`, or
#'   `"gray"`.
#' @param legend Legend position: `"bottom"`, `"top"`, `"left"`, `"right"`, or
#'   `"none"`.
#' @param plot_opts Nested list for advanced plot customization. See the dedicated
#'   Plot options section below. Values supplied here override defaults.
#' @param lang `"en"` or `"vi"` for generated labels/messages.
#' @param digit Number of decimals for ordinary estimates in Viewer tables.
#' @param p_digit Number of decimals for p-values in Viewer tables.
#' @param interpretation Logical; include a short interpretation section in the
#'   Viewer/console report. Default is `FALSE`; the interpretation text is still
#'   retained in `$interpretation` for programmatic use.
#' @param viewer_plot_format Self-contained Viewer image format: `"png"` or
#'   `"svg"`. Both are produced with base R graphics and need no extra package.
#' @param console Print a concise analysis summary to the console. Default `FALSE`.
#' @param show Open the complete HTML report in the RStudio Viewer (or browser) and
#'   display requested figures in the Plot pane. Default `TRUE`.
#'
#' @section Rolling-window interpretation:
#' With `window = 5`, `window_type = "rolling"`, `window_align = "right"`, and
#' `window_step = 1`, the first displayed point summarizes cases 1-5, the next
#' summarizes 2-6, then 3-7, and so on. Therefore the point at case 100 represents
#' performance in the most recent five cases (96-100). A newly observed case 101
#' updates the curve to cases 97-101. This is different from non-overlapping block
#' summaries and from a cumulative mean.
#'
#' @section Statistical inference versus visual smoothing:
#' Overlapping rolling windows share observations and are therefore correlated.
#' `tablearn()` can display rolling windows for a stable curve while fitting the
#' change-point model and CUSUM on the original case sequence. The recommended
#' default is `break_on = "raw"`; CUSUM, LC-CUSUM and RA-CUSUM always operate on
#' individual sequential cases in this implementation.
#'
#' @section Advanced plot options:
#' `plot_opts` is a nested list. Every field is optional. Main groups are:
#'
#' * `raw`: `show`, `color`, `fill`, `shape`, `size`, `alpha`, `stroke`.
#' * `window`: `show`, `geom` (`"line"`, `"point"`, `"point_line"`), `color`,
#'   `fill`, `shape`, `size`, `alpha`, `line_color`, `line_width`, `line_type`.
#' * `window_interval`: `show`, `geom` (`"ribbon"` or `"errorbar"`), `color`,
#'   `fill`, `alpha`, `line_width`, `width`.
#' * `smooth`: `show`, `color`, `fill`, `line_width`, `line_type`, `alpha`,
#'   `ci_show`, `ci_alpha`.
#' * `ewma`: `show`, `color`, `line_width`, `line_type`, `alpha`.
#' * `breakpoint`: `show`, `color`, `line_width`, `line_type`, `alpha`, `label`,
#'   `label_text`, `label_size`, `label_angle`, `label_hjust`, `label_vjust`.
#' * `proficiency`: independent proficiency-line controls: `show`, `color`,
#'   `line_width`, `line_type`, `alpha`, `label`, `label_text`, `label_size`,
#'   `label_angle`, `label_hjust`, `label_vjust`.
#' * `phase`: `show`, `fills`, `alpha`, `border_color`, `border_width`, `label`,
#'   `label_size`, `label_position`.
#' * `target`: `show`, `color`, `line_width`, `line_type`, `alpha`, `label`,
#'   `label_text`, `label_size`.
#' * `current`: `show`, `color`, `fill`, `shape`, `size`, `alpha`, `stroke`,
#'   `label`, `label_text`, `label_size`, `hjust`, `vjust`.
#' * `axes`: `x_limits`, `y_limits`, `x_breaks`, `y_breaks`, `x_expand`,
#'   `y_expand`, `y_percent`, `percent_accuracy`, `x_reverse`, `y_reverse`,
#'   `x_trans`, `y_trans`, `x_date_format`, `x_date_breaks`, `clip`.
#' * `legend`: `show`, `position`, `title`, `direction`.
#' * `facet`: `ncol`, `nrow`, `scales`.
#' * `operator`: `colors`, `shapes`, `line_types`; layout is selected by
#'   the high-level `operator_display` argument.
#' * `cusum`: CUSUM-figure controls including `colors`, `line_types`, `line_width`,
#'   `alpha`, point controls, zero-line controls, decision-limit controls, signal
#'   marker controls, title/subtitle/caption/axis labels, and CUSUM-specific axis
#'   limits/breaks.
#' * `theme`: `name`, `base_size`, `base_family`, `grid_major`, `grid_minor`,
#'   `panel_border`, `axis_line`, `plot_title_face`, `legend_key_size`.
#' * `text`: optional title/subtitle/caption/axis/legend text sizes.
#' * `margins`: `top`, `right`, `bottom`, `left` in points.
#' * `panel`: optional background/border customization.
#' * `annotation`: `show_n`, `show_window_label`.
#'
#' This layered design lets the raw observations remain visible while the rolling
#' curve, uncertainty, fitted smooth, target, breakpoint, estimated proficiency,
#' phases, and current performance are styled independently.
#'
#' @section Proficiency and phase classification:
#' `tablearn()` deliberately distinguishes a statistical/descriptive change point
#' from proficiency. A change point indicates a change in the learning trajectory;
#' proficiency is estimated from one or more operational criteria. With the default
#' combined rule, a supplied clinical target must be sustained for `target_hold`
#' consecutive displayed windows, and an enabled LC-CUSUM must cross its competency
#' decision boundary. A post-change plateau and reduced variability strengthen the
#' evidence. If no target or LC-CUSUM is supplied, a clear plateau can provide a
#' limited, data-driven proficiency estimate. Manual proficiency is also supported.
#'
#' Automatic phase classification uses these results rather than forcing every
#' dataset into three phases. When both an earlier learning change point and a later
#' proficiency case are found, phases are Learning -> Consolidation -> Proficiency.
#' If proficiency is not established, a post-change segment is labelled
#' Consolidation rather than Proficiency. `phase_breaks` always allows complete
#' manual control for study protocols with pre-specified phases.
#'
#' The evidence label (`Strong`, `Moderate`, `Limited`) is an R4VN rule-based
#' summary of concordant criteria, not a confidence probability and not a substitute
#' for a clinically defined competency standard.
#'
#' @section Viewer and dependency policy:
#' With `show = TRUE` (default), `tablearn()` opens one self-contained HTML report
#' containing the key publication-ready tables and every requested figure. The same
#' figures are also sent to the Plot pane. Standard analysis, HTML rendering, learning
#' curves, CUSUM, LC-CUSUM, and RA-CUSUM use only base/recommended R packages.
#' `ggplot2` is optional: when already installed, an advanced ggplot object is retained
#' in `$plots`; when it is absent, plotting still works through the base-R fallback.
#' `mgcv` is needed only when the user explicitly selects `smooth_method = "gam"`;
#' otherwise the default smoothing methods use base R.
#'
#' @return An object of class `r4vn_tablearn`. For one outcome it contains at least
#'   `raw`, `window`, `ewma`, `smooth`, `breakpoint`, `proficiency`,
#'   `proficiency_evidence`, `phases`, `phase_classification`, `current`,
#'   `current_status`, `cusum`, `lccusum`, `racusum`, `sensitivity`, `plots`,
#'   `settings`, `plot_options`, `tables`, and `interpretation`. With `show = TRUE`,
#'   the returned object also carries the generated self-contained Viewer HTML/file path.
#'   Multiple outcomes return class `r4vn_tablearn_multi` containing one analysis
#'   per outcome.
#'
#' @examples
#' # --------------------------------------------------------------------------
#' # Reproducible demonstration data used by the examples below
#' # --------------------------------------------------------------------------
#' set.seed(2026)
#' n <- 150
#' d <- data.frame(
#'   case = 1:n,
#'   date = as.Date("2025-01-01") + 0:(n - 1),
#'   surgeon = rep(c("A", "B", "C"), each = n / 3),
#'   complexity = rbinom(n, 1, 0.35),
#'   age = round(rnorm(n, 58, 12), 1)
#' )
#' d$time <- 115 - 48 * (1 - exp(-d$case / 30)) +
#'   9 * d$complexity + rnorm(n, 0, 9)
#' d$score <- 55 + 28 * (1 - exp(-d$case / 35)) + rnorm(n, 0, 5)
#' d$expected_risk <- plogis(-1.8 + 0.9 * d$complexity + 0.012 * (d$age - 58))
#' actual_risk <- plogis(qlogis(d$expected_risk) - 0.010 * d$case)
#' d$complication <- rbinom(n, 1, actual_risk)
#' d$errors <- rpois(n, pmax(0.15, 3.2 * exp(-d$case / 45)))
#'
#' # The full catalogue is interactive so R CMD check stays fast.
#' if (interactive()) {
#'
#' # 1. Simplest end-user command. In an interactive session this opens the
#' # complete Viewer report and sends the learning curve to the Plot pane.
#' if (interactive()) {
#'   m1 <- tablearn(time, data = d, order = case)
#' }
#'
#' # 2. Right-aligned rolling window: 1-5, 2-6, 3-7, ...
#' m2 <- tablearn(time, data = d, order = case, window = 5,
#'                show = FALSE, plot = FALSE)
#' head(m2$tables$Window_performance)
#' m2$tables$Current_performance
#'
#' # 3. Non-overlapping blocks: 1-10, 11-20, 21-30, ...
#' m3 <- tablearn(time, data = d, order = case, window = 10,
#'                window_type = "block", show = FALSE, plot = FALSE)
#' head(m3$window[, c(".start", ".end", ".value")])
#'
#' # 4. Cumulative learning curve: 1, 1-2, 1-3, ...
#' m4 <- tablearn(time, data = d, order = case,
#'                window_type = "cumulative", show = FALSE, plot = FALSE)
#'
#' # 5. Individual-case series without aggregation.
#' m5 <- tablearn(time, data = d, order = case,
#'                window_type = "none", show = FALSE, plot = FALSE)
#'
#' # 6. Median and IQR for a skewed continuous outcome.
#' m6 <- tablearn(time, data = d, order = case, window = 7,
#'                window_fun = "median", interval = "iqr",
#'                show = FALSE, plot = FALSE)
#'
#' # 7. Give recent cases greater weight within the rolling window.
#' m7 <- tablearn(time, data = d, order = case, window = 10,
#'                window_weight = "exponential", window_decay = 0.85,
#'                show = FALSE, plot = FALSE)
#'
#' # 8. Add EWMA to the ordinary rolling curve.
#' m8 <- tablearn(time, data = d, order = case, window = 10,
#'                ewma = TRUE, ewma_lambda = 0.20,
#'                show = FALSE, plot = FALSE)
#' tail(m8$ewma)
#'
#' # 9. Automatic piecewise change point. Inference uses raw cases by default.
#' m9 <- tablearn(time, data = d, order = case, window = 5,
#'                breakpoint = TRUE, break_on = "raw",
#'                show = FALSE, plot = FALSE)
#' m9$tables$Change_point
#'
#' # 10. Sustained clinical target: <= 70 for 3 consecutive windows.
#' m10 <- tablearn(time, data = d, order = case, window = 5,
#'                 better = "lower", target = 70, target_hold = 3,
#'                 proficiency_method = "target",
#'                 show = FALSE, plot = FALSE)
#' m10$tables$Proficiency
#' m10$tables$Proficiency_evidence
#'
#' # 11. Report the first qualifying window rather than the confirmation window.
#' m11 <- tablearn(time, data = d, order = case, window = 5,
#'                 better = "lower", target = 70, target_hold = 3,
#'                 proficiency_method = "target",
#'                 proficiency_case_rule = "first_stable",
#'                 show = FALSE, plot = FALSE)
#'
#' # 12. Manual proficiency and prespecified study phases.
#' m12 <- tablearn(time, data = d, order = case, window = 5,
#'                 proficiency_method = "manual", proficiency_case = 60,
#'                 phase_method = "manual", phase_breaks = c(25, 59),
#'                 phase_names = c("Learning", "Consolidation", "Proficiency"),
#'                 show = FALSE, plot = FALSE)
#' m12$tables$Phase_classification
#'
#' # 13. Binary outcome. Numeric 0/1 is recognized automatically; event = 1 is
#' # explicit and makes the scientific meaning clear.
#' m13 <- tablearn(complication, data = d, order = case, event = 1,
#'                 better = "lower", window = 20,
#'                 window_ci_method = "wilson",
#'                 show = FALSE, plot = FALSE)
#' m13$tables$Current_performance
#'
#' # 14. Count outcome.
#' m14 <- tablearn(errors, data = d, order = case, type = "count",
#'                 better = "lower", window = 10,
#'                 show = FALSE, plot = FALSE)
#'
#' # 15. Higher values can represent better performance.
#' m15 <- tablearn(score, data = d, order = case, better = "higher",
#'                 target = 80, target_hold = 3,
#'                 show = FALSE, plot = FALSE)
#'
#' # 16. Conventional deviation CUSUM for a continuous outcome.
#' m16 <- tablearn(time, data = d, order = case, target = 75,
#'                 better = "lower", cusum = TRUE,
#'                 plot_type = "auto", show = FALSE, plot = FALSE)
#' m16$tables$Sequential_monitoring
#'
#' # 17. Binary likelihood-ratio CUSUM.
#' m17 <- tablearn(complication, data = d, order = case, event = 1,
#'                 better = "lower", target = 0.10,
#'                 cusum = TRUE, cusum_method = "llr", cusum_alt = 0.20,
#'                 show = FALSE, plot = FALSE)
#'
#' # 18. LC-CUSUM: evidence that an acceptable failure rate has been reached.
#' m18 <- tablearn(complication, data = d, order = case, event = 1,
#'                 better = "lower", window = 20,
#'                 lccusum = TRUE, p_acceptable = 0.10,
#'                 p_unacceptable = 0.25,
#'                 proficiency_method = "lccusum",
#'                 show = FALSE, plot = FALSE)
#' m18$tables$Proficiency
#'
#' # 19. RA-CUSUM with externally supplied case-specific expected risks.
#' m19 <- tablearn(complication, data = d, order = case, event = 1,
#'                 better = "lower", window = 20,
#'                 racusum = TRUE, expected = expected_risk,
#'                 racusum_or = 2,
#'                 show = FALSE, plot = FALSE)
#'
#' # 20. RA-CUSUM can estimate expected risk from covariates using base glm().
#' # For prospective monitoring, an external/pre-specified risk model is preferred.
#' m20 <- suppressWarnings(tablearn(
#'   complication, data = d, order = case, event = 1,
#'   racusum = TRUE, adjust = vars(complexity, c.age), racusum_or = 2,
#'   show = FALSE, plot = FALSE
#' ))
#'
#' # 21. Separate learning curves by operator/surgeon.
#' m21 <- tablearn(time, data = d, order = case, operator = surgeon,
#'                 window = 8, operator_display = "facet",
#'                 show = FALSE, plot = FALSE)
#' m21$tables$Current_performance
#'
#' # 22. Overlay operators in the same graph.
#' m22 <- tablearn(time, data = d, order = case, operator = surgeon,
#'                 window = 8, operator_display = "overlay",
#'                 show = FALSE, plot = FALSE)
#'
#' # 23. Window-size sensitivity analysis.
#' m23 <- tablearn(time, data = d, order = case, window = 5,
#'                 sensitivity = TRUE,
#'                 window_sensitivity = c(3, 5, 10, 20),
#'                 show = FALSE, plot = FALSE)
#' m23$tables$Window_sensitivity
#'
#' # 24. Multiple outcomes in one command.
#' m24 <- tablearn(vars(time, complication), data = d, order = case,
#'                 event = 1, window = 10,
#'                 show = FALSE, plot = FALSE)
#' names(m24$outcomes)
#'
#' # 25. Use procedure date on the x-axis instead of consecutive case number.
#' m25 <- tablearn(time, data = d, order = date, x_axis = "order",
#'                 window = 7, show = FALSE, plot = FALSE)
#'
#' # 26. Interpretive prose is opt-in; default is FALSE.
#' m26 <- tablearn(time, data = d, order = case,
#'                 interpretation = TRUE, show = FALSE, plot = FALSE)
#' m26$interpretation
#'
#' # 27. Vietnamese generated interpretation/labels.
#' m27 <- tablearn(time, data = d, order = case, lang = "vi",
#'                 interpretation = TRUE, show = FALSE, plot = FALSE)
#'
#' # 28. Conservative auto-detection: a numeric variable with two values other
#' # than 0/1 remains continuous unless event/type explicitly says binary.
#' d2 <- data.frame(case = 1:20, value = c(rep(100, 10), rep(50, 10)))
#' m28 <- tablearn(value, data = d2, order = case,
#'                 show = FALSE, plot = FALSE)
#' m28$settings$type
#'
#' # 29. All publication-ready tables are directly accessible.
#' names(m10$tables)
#' summary(m10)$tables
#'
#' # 30. Plot methods work even when ggplot2 is not installed because tablearn()
#' # has a base-R plotting fallback. These also appear inside the Viewer report.
#' if (interactive()) {
#'   plot(m10, type = "learning")
#'   m30 <- tablearn(complication, data = d, order = case, event = 1,
#'                   lccusum = TRUE, p_acceptable = .10,
#'                   p_unacceptable = .25, plot_type = "both")
#'   plot(m30, type = "cusum")
#' }
#'
#' # 31. High-level publication styling. Advanced ggplot styling is used when
#' # ggplot2 is installed; the Viewer/base-R figure remains available otherwise.
#' if (interactive()) {
#'   m31 <- tablearn(
#'     time, data = d, order = case, window = 5, target = 70,
#'     plot_opts = list(
#'       raw = list(alpha = .15, size = 1.0),
#'       window = list(color = "#1F5A94", line_width = 1.2),
#'       smooth = list(color = "#B23A48", line_width = 1.4),
#'       proficiency = list(color = "#00796B"),
#'       phase = list(alpha = .08),
#'       current = list(fill = "#FFD166"),
#'       theme = list(base_size = 12, grid_minor = FALSE)
#'     )
#'   )
#' }
#'
#' } # end full interactive example catalogue
#'
#' @export
#' @md

tablearn <- function(
  outcome,
  data = NULL,
  order = NULL,
  operator = NULL,
  type = c("auto", "continuous", "binary", "count"),
  event = NULL,
  better = c("auto", "lower", "higher"),

  window = 5L,
  window_type = c("rolling", "block", "none", "cumulative"),
  window_align = c("right", "center", "left"),
  window_step = NULL,
  window_complete = TRUE,
  window_fun = c("auto", "mean", "median", "trimmed_mean", "proportion", "rate"),
  window_trim = 0.10,
  window_weight = c("equal", "linear", "exponential"),
  window_decay = 0.85,
  window_ci = TRUE,
  ci_level = 0.95,
  window_ci_method = c("auto", "t", "normal", "wilson", "exact"),
  interval = c("ci", "sd", "iqr", "range", "none"),

  ewma = FALSE,
  ewma_lambda = 0.20,
  ewma_init = c("first", "mean", "target"),

  smooth = TRUE,
  smooth_method = c("loess", "spline", "lm", "gam", "none"),
  smooth_on = c("window", "raw", "ewma"),
  smooth_span = 0.60,
  smooth_df = NULL,
  smooth_ci = TRUE,

  breakpoint = TRUE,
  break_on = c("raw", "window"),
  break_min_n = 8L,
  break_grid = NULL,
  break_boot = 0L,
  phase = TRUE,
  phase_breaks = NULL,
  phase_names = c("Learning", "Consolidation", "Proficiency"),
  phase_method = c("auto", "combined", "manual", "breakpoint", "proficiency"),

  proficiency = TRUE,
  proficiency_method = c("auto", "combined", "target", "lccusum", "breakpoint", "manual"),
  proficiency_case = NULL,
  proficiency_case_rule = c("confirmed", "first_stable"),
  target_hold = 3L,
  target_tolerance = 0,
  plateau = TRUE,
  plateau_ratio = 0.25,
  plateau_slope = NULL,
  stability = TRUE,
  stability_metric = c("auto", "sd", "iqr", "cv", "none"),
  stability_window = NULL,
  stability_ratio = 0.75,

  cusum = FALSE,
  target = NULL,
  cusum_method = c("deviation", "llr"),
  cusum_alt = NULL,
  cusum_reset = FALSE,
  cusum_limit = NULL,

  lccusum = FALSE,
  p_acceptable = NULL,
  p_unacceptable = NULL,
  alpha = 0.05,
  beta = 0.10,

  racusum = FALSE,
  expected = NULL,
  adjust = NULL,
  racusum_or = 2,
  racusum_limit = NULL,

  sensitivity = FALSE,
  window_sensitivity = NULL,

  plot = TRUE,
  plot_type = c("auto", "learning", "cusum", "both"),
  x_axis = c("case", "order"),
  operator_display = c("facet", "overlay"),
  show_raw = TRUE,
  show_window = TRUE,
  show_smooth = TRUE,
  show_ewma = TRUE,
  show_interval = TRUE,
  show_break = TRUE,
  show_proficiency = TRUE,
  show_phase = TRUE,
  show_target = TRUE,
  show_current = TRUE,
  title = NULL,
  subtitle = NULL,
  caption = NULL,
  xlab = NULL,
  ylab = NULL,
  theme = c("publication", "minimal", "classic", "bw", "gray"),
  legend = "bottom",
  plot_opts = NULL,
  lang = c("en", "vi"),
  digit = 2,
  p_digit = 3,
  interpretation = FALSE,
  viewer_plot_format = c("png", "svg"),
  console = FALSE,
  show = TRUE
) {
  # Resolve NSE-sensitive arguments immediately in the caller's environment.
  # This prevents promises such as `data = lc_test` from being re-evaluated later
  # inside do.call()/internal functions, where `lc_test` may no longer be visible.
  call <- match.call()
  caller_env <- parent.frame()

  plot_type <- match.arg(plot_type)
  if (identical(plot_type, "auto")) {
    plot_type <- if (isTRUE(cusum) || isTRUE(lccusum) || isTRUE(racusum)) "both" else "learning"
  }
  viewer_plot_format <- match.arg(viewer_plot_format)
  digit <- as.integer(digit)
  p_digit <- as.integer(p_digit)
  if (length(digit) != 1L || is.na(digit) || digit < 0L || digit > 10L) {
    .tablearn_stop("`digit` must be one integer from 0 to 10.")
  }
  if (length(p_digit) != 1L || is.na(p_digit) || p_digit < 1L || p_digit > 10L) {
    .tablearn_stop("`p_digit` must be one integer from 1 to 10.")
  }
  if (!is.logical(interpretation) || length(interpretation) != 1L || is.na(interpretation)) {
    .tablearn_stop("`interpretation` must be TRUE or FALSE.")
  }

  data_expr <- substitute(data)
  data_value <- if (missing(data) || identical(data_expr, quote(NULL))) {
    NULL
  } else {
    tryCatch(
      eval(data_expr, envir = caller_env),
      error = function(e) {
        .tablearn_stop(
          "Could not evaluate `data` in the calling environment: ",
          conditionMessage(e)
        )
      }
    )
  }
  data_value <- .tablearn_active_data(data_value)

  out_expr <- substitute(outcome)
  ord_expr <- substitute(order)
  op_expr <- substitute(operator)
  adj_expr <- substitute(adjust)
  exp_expr <- substitute(expected)

  outcomes <- .tablearn_extract_vars_from_expr(out_expr, data_value, caller_env, allow_null = FALSE)
  order_var <- .tablearn_extract_vars_from_expr(ord_expr, data_value, caller_env, allow_null = TRUE)
  operator_var <- .tablearn_extract_vars_from_expr(op_expr, data_value, caller_env, allow_null = TRUE)
  adjust_vars <- .tablearn_extract_vars_from_expr(adj_expr, data_value, caller_env, allow_null = TRUE)

  if (length(order_var) > 1L) .tablearn_stop("`order` must identify one variable.")
  if (length(operator_var) > 1L) .tablearn_stop("`operator` must identify one variable.")
  .tablearn_check_vars(outcomes, data_value, "outcome")
  .tablearn_check_vars(order_var, data_value, "order")
  .tablearn_check_vars(operator_var, data_value, "operator")
  .tablearn_check_vars(adjust_vars, data_value, "adjust")

  # Resolve expected risk now as well. A bare column name is retained as a column
  # reference; expressions such as `lc_test$risk` are converted to a concrete vector.
  expected_value <- NULL
  if (!identical(exp_expr, quote(NULL))) {
    if (is.symbol(exp_expr) && as.character(exp_expr) %in% names(data_value)) {
      expected_value <- as.character(exp_expr)
    } else {
      expected_value <- tryCatch(
        eval(exp_expr, envir = caller_env),
        error = function(e) {
          .tablearn_stop(
            "Could not evaluate `expected` in the calling environment: ",
            conditionMessage(e)
          )
        }
      )
    }
  }

  # Force ordinary arguments before entering the internal core. This makes the
  # internal implementation completely standard-evaluation based.
  resolved_args <- list(
    data = data_value,
    order_var = order_var,
    operator_var = operator_var,
    type = match.arg(type),
    event = event,
    better = match.arg(better),
    window = window,
    window_type = match.arg(window_type),
    window_align = match.arg(window_align),
    window_step = window_step,
    window_complete = window_complete,
    window_fun = match.arg(window_fun),
    window_trim = window_trim,
    window_weight = match.arg(window_weight),
    window_decay = window_decay,
    window_ci = window_ci,
    ci_level = ci_level,
    window_ci_method = match.arg(window_ci_method),
    interval = match.arg(interval),
    ewma = ewma,
    ewma_lambda = ewma_lambda,
    ewma_init = match.arg(ewma_init),
    smooth = smooth,
    smooth_method = match.arg(smooth_method),
    smooth_on = match.arg(smooth_on),
    smooth_span = smooth_span,
    smooth_df = smooth_df,
    smooth_ci = smooth_ci,
    breakpoint = breakpoint,
    break_on = match.arg(break_on),
    break_min_n = as.integer(break_min_n),
    break_grid = break_grid,
    break_boot = as.integer(break_boot),
    phase = phase,
    phase_breaks = phase_breaks,
    phase_names = phase_names,
    phase_method = match.arg(phase_method),
    proficiency = proficiency,
    proficiency_method = match.arg(proficiency_method),
    proficiency_case = proficiency_case,
    proficiency_case_rule = match.arg(proficiency_case_rule),
    target_hold = as.integer(target_hold),
    target_tolerance = target_tolerance,
    plateau = plateau,
    plateau_ratio = plateau_ratio,
    plateau_slope = plateau_slope,
    stability = stability,
    stability_metric = match.arg(stability_metric),
    stability_window = stability_window,
    stability_ratio = stability_ratio,
    cusum = cusum,
    target = target,
    cusum_method = match.arg(cusum_method),
    cusum_alt = cusum_alt,
    cusum_reset = cusum_reset,
    cusum_limit = cusum_limit,
    lccusum = lccusum,
    p_acceptable = p_acceptable,
    p_unacceptable = p_unacceptable,
    alpha = alpha,
    beta = beta,
    racusum = racusum,
    expected_value = expected_value,
    adjust_vars = adjust_vars,
    racusum_or = racusum_or,
    racusum_limit = racusum_limit,
    sensitivity = sensitivity,
    window_sensitivity = window_sensitivity,
    plot = plot,
    plot_type = plot_type,
    x_axis = match.arg(x_axis),
    operator_display = match.arg(operator_display),
    show_raw = show_raw,
    show_window = show_window,
    show_smooth = show_smooth,
    show_ewma = show_ewma,
    show_interval = show_interval,
    show_break = show_break,
    show_proficiency = show_proficiency,
    show_phase = show_phase,
    show_target = show_target,
    show_current = show_current,
    title = title,
    subtitle = subtitle,
    caption = caption,
    xlab = xlab,
    ylab = ylab,
    theme = match.arg(theme),
    legend = legend,
    plot_opts = plot_opts,
    lang = match.arg(lang)
  )

  run_one <- function(outcome_name) {
    # IMPORTANT: do not pass `call <- match.call()` through do.call().
    # A language object inside `args` is evaluated by do.call(quote = FALSE);
    # passing it would try to execute tablearn(...) again in `envir`, which is
    # exactly what caused "could not find function 'tablearn'" during tests.
    z <- do.call(
      .r4vn_tablearn_core,
      c(list(outcome_var = outcome_name), resolved_args),
      quote = FALSE
    )
    z$call <- call
    z$settings$plot <- isTRUE(plot)
    z$settings$digit <- digit
    z$settings$p_digit <- p_digit
    z$settings$interpretation <- isTRUE(interpretation)
    z$settings$viewer_plot_format <- viewer_plot_format
    z$settings$title <- title
    z$settings$subtitle <- subtitle
    z$settings$caption <- caption
    z$settings$plot_type <- plot_type
    z$tables <- .tablearn_make_tables(z, digit = digit, p_digit = p_digit)
    z
  }

  if (length(outcomes) > 1L) {
    ans <- lapply(outcomes, run_one)
    names(ans) <- outcomes
    out <- list(
      call = call, outcomes = ans, outcome_names = outcomes,
      settings = list(
        plot = isTRUE(plot), plot_type = plot_type, digit = digit,
        p_digit = p_digit, interpretation = isTRUE(interpretation),
        viewer_plot_format = viewer_plot_format, title = title
      )
    )
    class(out) <- c("r4vn_tablearn_multi", "list")

    out <- .r4vn_show(out, show = show, console = console, renderer = .r4vn_viewer_tablearn)
    if (isTRUE(show) && isTRUE(plot)) {
      invisible(lapply(ans, function(z) {
        if (plot_type %in% c("learning", "both")) plot(z, type = "learning")
        if (plot_type %in% c("cusum", "both") && .tablearn_has_cusum(z)) plot(z, type = "cusum")
        invisible(NULL)
      }))
    }
    return(invisible(out))
  }

  out <- run_one(outcomes[[1L]])
  out <- .r4vn_show(out, show = show, console = console, renderer = .r4vn_viewer_tablearn)
  if (isTRUE(show) && isTRUE(plot)) {
    if (plot_type %in% c("learning", "both")) plot(out, type = "learning")
    if (plot_type %in% c("cusum", "both") && .tablearn_has_cusum(out)) plot(out, type = "cusum")
  }
  invisible(out)
}

.r4vn_tablearn_core <- function(outcome_var, data, order_var, operator_var, type, event, better,
                           window, window_type, window_align, window_step, window_complete,
                           window_fun, window_trim, window_weight, window_decay, window_ci,
                           ci_level, window_ci_method, interval, ewma, ewma_lambda, ewma_init,
                           smooth, smooth_method, smooth_on, smooth_span, smooth_df, smooth_ci,
                           breakpoint, break_on, break_min_n, break_grid, break_boot, phase,
                           phase_breaks, phase_names, phase_method, proficiency,
                           proficiency_method, proficiency_case, proficiency_case_rule,
                           target_hold, target_tolerance, plateau, plateau_ratio, plateau_slope,
                           stability, stability_metric, stability_window, stability_ratio,
                           cusum, target, cusum_method, cusum_alt,
                           cusum_reset, cusum_limit, lccusum, p_acceptable, p_unacceptable,
                           alpha, beta, racusum, expected_value, adjust_vars, racusum_or,
                           racusum_limit, sensitivity, window_sensitivity, plot, plot_type,
                           x_axis, operator_display, show_raw, show_window, show_smooth,
                           show_ewma, show_interval, show_break, show_proficiency, show_phase, show_target,
                           show_current, title, subtitle, caption, xlab, ylab, theme, legend,
                           plot_opts, lang, call = NULL) {

  if (!is.numeric(window) || length(window) != 1L || window < 1) .tablearn_stop("`window` must be a positive integer.")
  window <- as.integer(window)
  if (is.null(window_step)) window_step <- if (window_type == "block") window else 1L
  window_step <- as.integer(window_step)
  if (window_step < 1L) .tablearn_stop("`window_step` must be >= 1.")
  if (!is.numeric(ci_level) || ci_level <= 0 || ci_level >= 1) .tablearn_stop("`ci_level` must lie in (0, 1).")
  if (length(target_hold) != 1L || is.na(target_hold) || target_hold < 1L) .tablearn_stop("`target_hold` must be one integer >= 1.")
  if (!is.null(stability_window)) {
    stability_window <- as.integer(stability_window)
    if (length(stability_window) != 1L || is.na(stability_window) || stability_window < 3L) .tablearn_stop("`stability_window` must be NULL or one integer >= 3.")
  }

  y0 <- data[[outcome_var]]
  outcome_type <- .tablearn_detect_type(y0, type, event)
  event_used <- event
  if (window_fun == "proportion" && outcome_type != "binary") {
    .tablearn_stop("`window_fun = 'proportion'` requires a binary outcome.")
  }
  if (outcome_type == "binary" && interval %in% c("sd", "iqr", "range")) {
    .tablearn_stop("For a binary outcome, use `interval = 'ci'` or `interval = 'none'`.")
  }
  if (outcome_type == "binary" && window_ci_method == "t") {
    .tablearn_stop("`window_ci_method = 't'` is not appropriate for binary proportions. Use 'wilson', 'exact', 'normal', or 'auto'.")
  }
  if (outcome_type != "binary" && window_ci_method %in% c("wilson", "exact")) {
    .tablearn_stop("Wilson/exact intervals are for binary proportions. Use 't', 'normal', or 'auto' for continuous outcomes.")
  }
  if (outcome_type == "binary") {
    bb <- .tablearn_binary(y0, event)
    y_analysis <- bb$value
    event_used <- bb$event
  } else {
    if (!is.numeric(y0) && !is.integer(y0)) .tablearn_stop("Continuous/count outcomes must be numeric.")
    y_analysis <- as.numeric(y0)
    if (outcome_type == "count" && any(y_analysis < 0 | abs(y_analysis - round(y_analysis)) > sqrt(.Machine$double.eps), na.rm = TRUE)) {
      .tablearn_stop("`type = 'count'` requires non-negative integer outcomes.")
    }
  }

  if (better == "auto") better <- "lower"
  if (window_fun %in% c("median", "trimmed_mean") && window_weight != "equal") {
    .tablearn_warn("`window_weight` is ignored by median/trimmed-mean summaries; equal weighting is used for that summary statistic.")
  }
  if (isTRUE(breakpoint) && break_on == "window" && window_type == "rolling" && window_step < window) {
    .tablearn_warn("`break_on = 'window'` with overlapping rolling windows treats correlated summaries as model observations. For inference, `break_on = 'raw'` is recommended.")
  }

  ord <- if (is.null(order_var)) seq_len(nrow(data)) else data[[order_var]]
  if (x_axis == "order" && !(is.numeric(ord) || inherits(ord, "Date"))) {
    .tablearn_warn("`x_axis = 'order'` currently supports numeric or Date ordering variables; falling back to case number.")
    x_axis <- "case"
  }
  op <- if (is.null(operator_var)) rep("Overall", nrow(data)) else as.character(data[[operator_var]])
  op[is.na(op)] <- "(Missing operator)"

  work <- data.frame(.row = seq_len(nrow(data)), .operator = op, stringsAsFactors = FALSE)
  work$.order <- ord
  work$.value <- y_analysis
  keep <- !is.na(work$.value) & !is.na(work$.order)
  work <- work[keep, , drop = FALSE]
  if (!nrow(work)) .tablearn_stop("No complete outcome/order observations are available.")

  operator_levels <- unique(work$.operator)
  byop <- vector("list", length(operator_levels)); names(byop) <- operator_levels

  for (oo in operator_levels) {
    w <- work[work$.operator == oo, , drop = FALSE]
    w <- w[order(w$.order, w$.row), , drop = FALSE]
    w$.case <- seq_len(nrow(w))

    win <- .tablearn_compute_windows(
      y = w$.value, x_order = w$.order, outcome_type = outcome_type,
      size = window, type = window_type, align = window_align, step = window_step,
      complete = window_complete, fun = window_fun, trim = window_trim,
      weight = window_weight, decay = window_decay, ci = window_ci,
      ci_level = ci_level, ci_method = window_ci_method, interval = interval
    )
    if (nrow(win)) win$.operator <- oo

    ew <- if (isTRUE(ewma)) .tablearn_ewma(w$.value, ewma_lambda, ewma_init, target) else data.frame()
    if (nrow(ew)) {
      ew$.order <- w$.order[ew$.case]
      ew$.operator <- oo
    }

    sm <- data.frame()
    if (isTRUE(smooth) && smooth_method != "none") {
      if (smooth_on == "raw") { sx <- w$.case; sy <- w$.value }
      else if (smooth_on == "ewma" && nrow(ew)) { sx <- ew$.case; sy <- ew$.value }
      else if (nrow(win)) { sx <- win$.x_case; sy <- win$.value }
      else { sx <- w$.case; sy <- w$.value }
      sm <- .tablearn_smooth(sx, sy, smooth_method, smooth_span, smooth_df, smooth_ci, ci_level, outcome_type)
      if (nrow(sm)) sm$.operator <- oo
    }

    bp <- NULL
    if (isTRUE(breakpoint)) {
      if (break_on == "window" && nrow(win)) { bx <- win$.x_case; by <- win$.value }
      else { bx <- w$.case; by <- w$.value }
      bp <- .tablearn_breakpoint(bx, by, outcome_type, break_min_n, break_grid, break_boot, ci_level)
      if (!is.null(bp)) bp$operator <- oo
    }

    cs <- if (isTRUE(cusum)) .tablearn_cusum(w$.value, outcome_type, target, cusum_alt, better, cusum_method, cusum_reset, cusum_limit) else data.frame()
    if (nrow(cs)) cs$.operator <- oo

    lc <- if (isTRUE(lccusum)) {
      if (outcome_type != "binary") .tablearn_stop("LC-CUSUM requires a binary outcome.")
      .tablearn_lccusum(w$.value, p_acceptable, p_unacceptable, alpha, beta)
    } else data.frame()
    if (nrow(lc)) lc$.operator <- oo

    rac <- data.frame()
    if (isTRUE(racusum)) {
      if (outcome_type != "binary") .tablearn_stop("RA-CUSUM requires a binary outcome.")
      original_rows <- w$.row
      pexp <- NULL
      if (is.character(expected_value) && length(expected_value) == 1L && expected_value %in% names(data)) {
        pexp <- as.numeric(data[[expected_value]][original_rows])
      } else if (is.numeric(expected_value)) {
        if (length(expected_value) == nrow(data)) pexp <- expected_value[original_rows]
        else if (length(expected_value) == nrow(w)) pexp <- expected_value
        else .tablearn_stop("Numeric `expected` must have length nrow(data) or the number of analyzed cases.")
      } else if (length(adjust_vars)) {
        .tablearn_warn("RA-CUSUM expected risks are being estimated from the same analyzed dataset. For prospective monitoring, externally validated or pre-specified expected risks are preferable.")
        dd <- data[original_rows, c(outcome_var, adjust_vars), drop = FALSE]
        dd$.event <- w$.value
        frm <- stats::as.formula(paste(".event ~", paste(sprintf("`%s`", adjust_vars), collapse = " + ")))
        fitrisk <- stats::glm(frm, family = stats::binomial(), data = dd)
        pexp <- stats::predict(fitrisk, newdata = dd, type = "response")
      } else {
        .tablearn_stop("RA-CUSUM requires `expected` risks or `adjust` covariates.")
      }
      rac <- .tablearn_racusum(w$.value, pexp, racusum_or, racusum_limit)
      rac$.operator <- oo
    }

    prof_obj <- .tablearn_proficiency_engine(
      raw = w, win = win, bp = bp, lccusum = lc, racusum = rac,
      outcome_type = outcome_type, better = better, target = target,
      proficiency = proficiency, proficiency_method = proficiency_method,
      proficiency_case = proficiency_case, proficiency_case_rule = proficiency_case_rule,
      target_hold = target_hold, target_tolerance = target_tolerance,
      plateau = plateau, plateau_ratio = plateau_ratio, plateau_slope = plateau_slope,
      stability = stability, stability_metric = stability_metric,
      stability_window = stability_window, stability_ratio = stability_ratio,
      window_n = window, operator = oo
    )
    prof <- prof_obj$summary
    prof_ev <- prof_obj$evidence
    current_status <- prof_obj$current_status

    brs <- phase_breaks
    phase_names_used <- phase_names
    pm <- if (phase_method == "auto") "combined" else phase_method
    if (isTRUE(phase)) {
      if (!is.null(phase_breaks)) {
        brs <- phase_breaks
        phase_names_used <- phase_names
      } else if (pm == "manual") {
        .tablearn_stop("`phase_method = 'manual'` requires `phase_breaks`.")
      } else if (pm == "breakpoint") {
        if (!is.null(bp) && isTRUE(bp$detected)) {
          brs <- bp$breakpoint
          phase_names_used <- c("Learning", "Consolidation")
        } else {
          brs <- NULL
        }
      } else if (pm == "proficiency") {
        pc <- prof$proficiency_case[1L]
        if (isTRUE(prof$achieved[1L]) && is.finite(pc) && pc > min(w$.case)) {
          brs <- pc - 1
          phase_names_used <- c("Learning", "Proficiency")
        } else {
          brs <- NULL
        }
      } else {
        brs <- prof_obj$phase_breaks
        auto_names <- prof_obj$phase_names
        if (length(auto_names) == 3L && length(phase_names) >= 3L) {
          phase_names_used <- phase_names[1:3]
        } else if (length(auto_names) == 2L && length(phase_names) >= 2L) {
          if (isTRUE(prof$achieved[1L]) && length(phase_names) >= 3L) {
            phase_names_used <- c(phase_names[1L], phase_names[3L])
          } else {
            phase_names_used <- phase_names[1:2]
          }
        } else {
          phase_names_used <- auto_names
        }
      }
    }
    ph <- if (isTRUE(phase) && !is.null(brs) && length(brs)) {
      .tablearn_phase_table(w, brs, outcome_type, phase_names_used, ci_level)
    } else data.frame()
    if (nrow(ph)) ph$operator <- oo

    cur <- if (nrow(win)) {
      last <- win[nrow(win), , drop = FALSE]
      prev <- if (nrow(win) >= 2L) win$.value[nrow(win) - 1L] else NA_real_
      data.frame(
        operator = oo, case = last$.end, window_start = last$.start, window_end = last$.end,
        n = last$.n, events = last$.events, estimate = last$.value,
        lower = last$.lower, upper = last$.upper, previous = prev,
        change = last$.value - prev,
        target = if (is.null(target)) NA_real_ else target,
        target_met = if (is.null(target)) NA else if (better == "higher") last$.value >= target else last$.value <= target,
        stringsAsFactors = FALSE
      )
    } else data.frame()

    sens <- data.frame()
    if (isTRUE(sensitivity)) {
      ws <- window_sensitivity
      if (is.null(ws)) {
        nn <- nrow(w)
        ws <- unique(pmax(2L, pmin(nn, c(3L, 5L, 10L, 20L))))
      }
      sr <- lapply(ws, function(k) {
        wi <- .tablearn_compute_windows(w$.value, w$.order, outcome_type, as.integer(k), "rolling", "right", 1L, TRUE,
                                        window_fun, window_trim, window_weight, window_decay, FALSE, ci_level,
                                        window_ci_method, "none")
        bpw <- if (nrow(wi) >= 2L * break_min_n + 2L) .tablearn_breakpoint(wi$.x_case, wi$.value, outcome_type, break_min_n, break_grid, 0L, ci_level) else NULL
        data.frame(window = k, points = nrow(wi), current = if (nrow(wi)) tail(wi$.value, 1L) else NA_real_,
                   breakpoint = if (is.null(bpw)) NA_real_ else bpw$breakpoint,
                   detected = if (is.null(bpw)) NA else bpw$detected)
      })
      sens <- do.call(rbind, sr); sens$operator <- oo
    }

    byop[[oo]] <- list(raw = w, window = win, ewma = ew, smooth = sm, breakpoint = bp,
                       proficiency = prof, proficiency_evidence = prof_ev, phases = ph,
                       phase_classification = ph, current_status = current_status,
                       cusum = cs, lccusum = lc, racusum = rac,
                       current = cur, sensitivity = sens)
  }

  bind <- function(name) {
    xs <- lapply(byop, `[[`, name)
    xs <- xs[vapply(xs, function(z) is.data.frame(z) && nrow(z) > 0L, logical(1))]
    if (!length(xs)) return(data.frame())
    do.call(rbind, xs)
  }
  raw_all <- bind("raw"); win_all <- bind("window"); ew_all <- bind("ewma"); sm_all <- bind("smooth")
  prof_all <- bind("proficiency"); prof_ev_all <- bind("proficiency_evidence")
  ph_all <- bind("phases"); cs_all <- bind("cusum"); lc_all <- bind("lccusum"); ra_all <- bind("racusum")
  cur_all <- bind("current"); current_status_all <- bind("current_status"); sens_all <- bind("sensitivity")
  bp_all <- lapply(byop, `[[`, "breakpoint")

  interpretation <- .tablearn_interpret(outcome_var, outcome_type, better, byop, window, window_type, lang)

  settings <- list(
    outcome = outcome_var, order = order_var, operator = operator_var, type = outcome_type,
    event = event_used, better = better, window = window, window_type = window_type,
    window_align = window_align, window_step = window_step, window_fun = window_fun,
    window_weight = window_weight, ci_level = ci_level, target = target,
    ewma = isTRUE(ewma), smooth = isTRUE(smooth), smooth_method = smooth_method, smooth_on = smooth_on,
    breakpoint = isTRUE(breakpoint), break_on = break_on, phase = isTRUE(phase),
    phase_method = phase_method, proficiency = proficiency,
    proficiency_method = proficiency_method, proficiency_case_rule = proficiency_case_rule,
    target_hold = target_hold, target_tolerance = target_tolerance, plateau = plateau,
    plateau_ratio = plateau_ratio, plateau_slope = plateau_slope, stability = stability,
    stability_metric = stability_metric, stability_window = stability_window,
    stability_ratio = stability_ratio,
    cusum = isTRUE(cusum), lccusum = isTRUE(lccusum), racusum = isTRUE(racusum),
    sensitivity = isTRUE(sensitivity), cusum_limit = cusum_limit, racusum_limit = racusum_limit,
    p_acceptable = p_acceptable, p_unacceptable = p_unacceptable,
    operator_display = operator_display, x_axis = x_axis, lang = lang
  )

  po <- .tablearn_recursive_modify(.tablearn_plot_defaults(), plot_opts)
  po$raw$show <- isTRUE(show_raw) && isTRUE(po$raw$show)
  po$window$show <- isTRUE(show_window) && isTRUE(po$window$show)
  po$smooth$show <- isTRUE(show_smooth) && isTRUE(po$smooth$show)
  po$ewma$show <- isTRUE(show_ewma) && isTRUE(po$ewma$show)
  po$window_interval$show <- isTRUE(show_interval) && isTRUE(po$window_interval$show)
  po$breakpoint$show <- isTRUE(show_break) && isTRUE(po$breakpoint$show)
  po$proficiency$show <- isTRUE(show_proficiency) && isTRUE(po$proficiency$show)
  po$phase$show <- isTRUE(show_phase) && isTRUE(po$phase$show)
  po$target$show <- isTRUE(show_target) && isTRUE(po$target$show)
  po$current$show <- isTRUE(show_current) && isTRUE(po$current$show)
  po$operator$display <- operator_display
  po$theme$name <- theme
  po$legend$position <- legend

  out <- list(call = call, raw = raw_all, window = win_all, ewma = ew_all, smooth = sm_all,
              breakpoint = bp_all, proficiency = prof_all, proficiency_evidence = prof_ev_all,
              phases = ph_all, phase_classification = ph_all, current = cur_all,
              current_status = current_status_all, cusum = cs_all,
              lccusum = lc_all, racusum = ra_all, sensitivity = sens_all,
              by_operator = byop, settings = settings, plot_options = po,
              interpretation = interpretation, plots = list(learning = NULL, cusum = NULL))
  class(out) <- c("r4vn_tablearn", "list")

  if (isTRUE(plot) && requireNamespace("ggplot2", quietly = TRUE)) {
    if (plot_type %in% c("learning", "both")) {
      out$plots$learning <- .tablearn_plot_learning(out, po, title, subtitle, caption, xlab, ylab)
    }
    if (plot_type %in% c("cusum", "both")) {
      out$plots$cusum <- .tablearn_plot_cusum(out, po)
    }
  }

  out
}

.tablearn_interpret <- function(outcome, type, better, byop, window, window_type, lang) {
  one <- function(z, nm) {
    bp <- z$breakpoint
    cur <- z$current
    pr <- z$proficiency
    if (lang == "vi") {
      s <- paste0("Outcome `", outcome, "` cho ", nm, ": \u0111\u01b0\u1eddng cong \u0111\u01b0\u1ee3c m\u00f4 t\u1ea3 b\u1eb1ng ", window_type,
                  " window n=", window, ".")
      if (!is.null(bp)) {
        if (isTRUE(bp$detected)) s <- paste0(s, " M\u00f4 h\u00ecnh piecewise ghi nh\u1eadn \u0111i\u1ec3m thay \u0111\u1ed5i kho\u1ea3ng ca ", round(bp$breakpoint, 1), ".")
        else s <- paste0(s, " Ch\u01b0a c\u00f3 b\u1eb1ng ch\u1ee9ng BIC r\u00f5 cho m\u1ed9t \u0111i\u1ec3m thay \u0111\u1ed5i duy nh\u1ea5t.")
      }
      if (is.data.frame(pr) && nrow(pr)) {
        if (isTRUE(pr$achieved[1])) {
          s <- paste0(s, " Proficiency \u0111\u01b0\u1ee3c \u01b0\u1edbc t\u00ednh \u0111\u1ea1t t\u1ea1i kho\u1ea3ng ca ", round(pr$proficiency_case[1], 1),
                      " (m\u1ee9c b\u1eb1ng ch\u1ee9ng theo quy t\u1eafc R4VN: ", pr$evidence[1], ").")
        } else if (identical(pr$status[1], "Not achieved")) {
          s <- paste0(s, " C\u00e1c ti\u00eau ch\u00ed proficiency \u0111\u00e3 ch\u1ecdn ch\u01b0a \u0111\u01b0\u1ee3c \u0111\u00e1p \u1ee9ng.")
        } else if (!identical(pr$status[1], "Not assessed")) {
          s <- paste0(s, " D\u1eef li\u1ec7u hi\u1ec7n ch\u01b0a \u0111\u1ee7 \u0111\u1ec3 x\u00e1c l\u1eadp proficiency theo quy t\u1eafc \u0111\u00e3 ch\u1ecdn.")
        }
      }
      if (nrow(cur)) s <- paste0(s, " Performance hi\u1ec7n t\u1ea1i c\u1ee7a window g\u1ea7n nh\u1ea5t = ", signif(cur$estimate[1], 4), ".")
    } else {
      s <- paste0("Outcome `", outcome, "` for ", nm, " was summarized using a ", window_type,
                  " window of ", window, " cases.")
      if (!is.null(bp)) {
        if (isTRUE(bp$detected)) s <- paste0(s, " Piecewise regression identified a change point near case ", round(bp$breakpoint, 1), ".")
        else s <- paste0(s, " BIC did not clearly support a single change point.")
      }
      if (is.data.frame(pr) && nrow(pr)) {
        if (isTRUE(pr$achieved[1])) {
          s <- paste0(s, " Proficiency was estimated at approximately case ", round(pr$proficiency_case[1], 1),
                      " (R4VN rule-based evidence: ", pr$evidence[1], ").")
        } else if (identical(pr$status[1], "Not achieved")) {
          s <- paste0(s, " The selected proficiency criteria were not met.")
        } else if (!identical(pr$status[1], "Not assessed")) {
          s <- paste0(s, " Available data were insufficient to establish proficiency under the selected rule.")
        }
      }
      if (nrow(cur)) s <- paste0(s, " Current-window performance = ", signif(cur$estimate[1], 4), ".")
    }
    s
  }
  paste(vapply(names(byop), function(nm) one(byop[[nm]], nm), character(1)), collapse = "\n")
}

.tablearn_theme <- function(po) {
  if (!requireNamespace("ggplot2", quietly = TRUE)) return(NULL)
  nm <- .tablearn_or(po$theme$name, "publication")
  base_size <- .tablearn_or(po$theme$base_size, 11)
  base_family <- .tablearn_or(po$theme$base_family, "")
  th <- switch(nm,
    minimal = ggplot2::theme_minimal(base_size = base_size, base_family = base_family),
    classic = ggplot2::theme_classic(base_size = base_size, base_family = base_family),
    bw = ggplot2::theme_bw(base_size = base_size, base_family = base_family),
    gray = ggplot2::theme_gray(base_size = base_size, base_family = base_family),
    ggplot2::theme_bw(base_size = base_size, base_family = base_family)
  )
  if (!isTRUE(po$theme$grid_major)) th <- th + ggplot2::theme(panel.grid.major = ggplot2::element_blank())
  if (!isTRUE(po$theme$grid_minor)) th <- th + ggplot2::theme(panel.grid.minor = ggplot2::element_blank())
  if (isTRUE(po$theme$panel_border)) th <- th + ggplot2::theme(panel.border = ggplot2::element_rect(fill = NA, colour = po$panel$border_color, linewidth = po$panel$border_width))
  if (isTRUE(po$theme$axis_line)) th <- th + ggplot2::theme(axis.line = ggplot2::element_line())
  if (!is.null(po$panel$background)) th <- th + ggplot2::theme(panel.background = ggplot2::element_rect(fill = po$panel$background, colour = NA))
  th <- th + ggplot2::theme(
    plot.title = ggplot2::element_text(face = .tablearn_or(po$theme$plot_title_face, "bold")),
    legend.position = if (!isTRUE(po$legend$show) || identical(po$legend$position, "none")) "none" else po$legend$position,
    legend.direction = .tablearn_or(po$legend$direction, "horizontal"),
    legend.key.size = grid::unit(.tablearn_or(po$theme$legend_key_size, 0.9), "lines"),
    plot.margin = ggplot2::margin(po$margins$top, po$margins$right, po$margins$bottom, po$margins$left, unit = "pt")
  )
  ts <- po$text
  if (!is.null(ts$title_size)) th <- th + ggplot2::theme(plot.title = ggplot2::element_text(size = ts$title_size, face = .tablearn_or(po$theme$plot_title_face, "bold")))
  if (!is.null(ts$subtitle_size)) th <- th + ggplot2::theme(plot.subtitle = ggplot2::element_text(size = ts$subtitle_size))
  if (!is.null(ts$caption_size)) th <- th + ggplot2::theme(plot.caption = ggplot2::element_text(size = ts$caption_size))
  if (!is.null(ts$axis_title_size)) th <- th + ggplot2::theme(axis.title = ggplot2::element_text(size = ts$axis_title_size))
  if (!is.null(ts$axis_text_size)) th <- th + ggplot2::theme(axis.text = ggplot2::element_text(size = ts$axis_text_size))
  if (!is.null(ts$legend_text_size)) th <- th + ggplot2::theme(legend.text = ggplot2::element_text(size = ts$legend_text_size))
  th
}

.tablearn_map_smooth_order <- function(sm, raw) {
  if (!nrow(sm)) return(sm)
  pieces <- lapply(split(sm, sm$.operator), function(z) {
    rr <- raw[raw$.operator == z$.operator[1L], , drop = FALSE]
    if (!nrow(rr)) { z$.x <- z$.case; return(z) }
    if (inherits(rr$.order, "Date")) {
      xx <- stats::approx(rr$.case, as.numeric(rr$.order), xout = z$.case, rule = 2)$y
      z$.x <- as.Date(xx, origin = "1970-01-01")
    } else if (is.numeric(rr$.order)) {
      z$.x <- stats::approx(rr$.case, rr$.order, xout = z$.case, rule = 2)$y
    } else {
      z$.x <- z$.case
    }
    z
  })
  do.call(rbind, pieces)
}

.tablearn_plot_learning <- function(x, po, title, subtitle, caption, xlab, ylab) {
  s <- x$settings
  raw <- x$raw; win <- x$window; sm <- x$smooth; ew <- x$ewma; cur <- x$current; prof <- x$proficiency
  multi_op <- length(unique(raw$.operator)) > 1L
  use_order <- s$x_axis == "order"

  if (use_order) {
    raw$.x <- raw$.order
    if (nrow(win)) win$.x <- win$.x_order
    if (nrow(ew)) ew$.x <- ew$.order
    if (nrow(sm)) sm <- .tablearn_map_smooth_order(sm, raw)
  } else {
    raw$.x <- raw$.case
    if (nrow(win)) win$.x <- win$.x_case
    if (nrow(ew)) ew$.x <- ew$.case
    if (nrow(sm)) sm$.x <- sm$.case
  }

  p <- ggplot2::ggplot()

  # Phase shading. With multiple operators, phase-specific rectangles work cleanly
  # in facets. They are intentionally suppressed in overlay mode because different
  # operators may have different breakpoints and overlapping backgrounds become
  # visually ambiguous.
  if (isTRUE(po$phase$show) && nrow(x$phases) && !use_order && !(multi_op && po$operator$display == "overlay")) {
    ph <- x$phases
    ph$.xmin <- ph$start - 0.5
    ph$.xmax <- ph$end + 0.5
    ph$.xmid <- (ph$start + ph$end) / 2
    ph$.phase_fill <- rep(po$phase$fills, length.out = nrow(ph))
    p <- p + ggplot2::geom_rect(
      data = ph,
      ggplot2::aes(xmin = .xmin, xmax = .xmax, ymin = -Inf, ymax = Inf, fill = .phase_fill),
      inherit.aes = FALSE,
      alpha = po$phase$alpha,
      colour = po$phase$border_color,
      linewidth = po$phase$border_width
    ) + ggplot2::scale_fill_identity(guide = "none")
    if (isTRUE(po$phase$label)) {
      phase_bottom <- identical(po$phase$label_position, "bottom")
      p <- p + ggplot2::geom_text(
        data = ph,
        ggplot2::aes(x = .xmid, y = if (phase_bottom) -Inf else Inf, label = phase),
        inherit.aes = FALSE,
        size = po$phase$label_size,
        vjust = if (phase_bottom) -0.4 else 1.4,
        colour = "#333333"
      )
    }
  }

  # Raw observations.
  if (isTRUE(po$raw$show) && nrow(raw)) {
    if (multi_op && po$operator$display == "overlay") {
      p <- p + ggplot2::geom_point(
        data = raw,
        ggplot2::aes(x = .x, y = .value, color = .operator, fill = .operator, shape = .operator, group = .operator),
        size = po$raw$size, alpha = po$raw$alpha, stroke = po$raw$stroke
      )
    } else {
      p <- p + ggplot2::geom_point(
        data = raw, ggplot2::aes(x = .x, y = .value),
        color = po$raw$color, fill = po$raw$fill, shape = po$raw$shape,
        size = po$raw$size, alpha = po$raw$alpha, stroke = po$raw$stroke
      )
    }
  }

  # Window interval.
  if (isTRUE(po$window_interval$show) && nrow(win) && any(is.finite(win$.lower)) && any(is.finite(win$.upper))) {
    if (identical(po$window_interval$geom, "errorbar")) {
      if (multi_op && po$operator$display == "overlay") {
        p <- p + ggplot2::geom_errorbar(
          data = win,
          ggplot2::aes(x = .x, ymin = .lower, ymax = .upper, color = .operator, group = .operator),
          alpha = po$window_interval$alpha, linewidth = po$window_interval$line_width,
          width = po$window_interval$width
        )
      } else {
        p <- p + ggplot2::geom_errorbar(
          data = win, ggplot2::aes(x = .x, ymin = .lower, ymax = .upper),
          color = po$window_interval$color, alpha = po$window_interval$alpha,
          linewidth = po$window_interval$line_width, width = po$window_interval$width
        )
      }
    } else {
      if (multi_op && po$operator$display == "overlay") {
        p <- p + ggplot2::geom_ribbon(
          data = win,
          ggplot2::aes(x = .x, ymin = .lower, ymax = .upper, fill = .operator, group = .operator),
          alpha = po$window_interval$alpha, colour = NA
        )
      } else {
        p <- p + ggplot2::geom_ribbon(
          data = win, ggplot2::aes(x = .x, ymin = .lower, ymax = .upper, group = .operator),
          fill = po$window_interval$fill, alpha = po$window_interval$alpha, colour = NA
        )
      }
    }
  }

  # Rolling/block/cumulative window line and points.
  if (isTRUE(po$window$show) && nrow(win)) {
    geom <- .tablearn_or(po$window$geom, "point_line")
    if (geom %in% c("line", "point_line")) {
      if (multi_op && po$operator$display == "overlay") {
        p <- p + ggplot2::geom_line(
          data = win,
          ggplot2::aes(x = .x, y = .value, color = .operator, linetype = .operator, group = .operator),
          linewidth = po$window$line_width, alpha = po$window$alpha
        )
      } else {
        p <- p + ggplot2::geom_line(
          data = win, ggplot2::aes(x = .x, y = .value, group = .operator),
          color = po$window$line_color, linewidth = po$window$line_width,
          linetype = po$window$line_type, alpha = po$window$alpha
        )
      }
    }
    if (geom %in% c("point", "point_line")) {
      if (multi_op && po$operator$display == "overlay") {
        p <- p + ggplot2::geom_point(
          data = win,
          ggplot2::aes(x = .x, y = .value, color = .operator, fill = .operator, shape = .operator, group = .operator),
          size = po$window$size, alpha = po$window$alpha, stroke = 0.5
        )
      } else {
        p <- p + ggplot2::geom_point(
          data = win, ggplot2::aes(x = .x, y = .value),
          color = po$window$color, fill = po$window$fill, shape = po$window$shape,
          size = po$window$size, alpha = po$window$alpha, stroke = 0.5
        )
      }
    }
    if (isTRUE(po$annotation$show_n)) {
      p <- p + ggplot2::geom_text(
        data = win, ggplot2::aes(x = .x, y = .value, label = paste0("n=", .n)),
        size = 2.8, vjust = -0.9, check_overlap = TRUE
      )
    } else if (isTRUE(po$annotation$show_window_label)) {
      p <- p + ggplot2::geom_text(
        data = win, ggplot2::aes(x = .x, y = .value, label = .window_label),
        size = 2.7, vjust = -0.9, check_overlap = TRUE
      )
    }
  }

  # Smoothed trend.
  if (isTRUE(po$smooth$show) && nrow(sm)) {
    if (isTRUE(po$smooth$ci_show) && any(is.finite(sm$.lower)) && any(is.finite(sm$.upper))) {
      if (multi_op && po$operator$display == "overlay") {
        p <- p + ggplot2::geom_ribbon(
          data = sm,
          ggplot2::aes(x = .x, ymin = .lower, ymax = .upper, fill = .operator, group = .operator),
          alpha = po$smooth$ci_alpha, colour = NA
        )
      } else {
        p <- p + ggplot2::geom_ribbon(
          data = sm, ggplot2::aes(x = .x, ymin = .lower, ymax = .upper, group = .operator),
          fill = po$smooth$fill, alpha = po$smooth$ci_alpha, colour = NA
        )
      }
    }
    if (multi_op && po$operator$display == "overlay") {
      p <- p + ggplot2::geom_line(
        data = sm, ggplot2::aes(x = .x, y = .value, color = .operator, group = .operator),
        linewidth = po$smooth$line_width, linetype = po$smooth$line_type, alpha = po$smooth$alpha
      )
    } else {
      p <- p + ggplot2::geom_line(
        data = sm, ggplot2::aes(x = .x, y = .value, group = .operator),
        color = po$smooth$color, linewidth = po$smooth$line_width,
        linetype = po$smooth$line_type, alpha = po$smooth$alpha
      )
    }
  }

  # EWMA.
  if (isTRUE(po$ewma$show) && nrow(ew)) {
    if (multi_op && po$operator$display == "overlay") {
      p <- p + ggplot2::geom_line(
        data = ew, ggplot2::aes(x = .x, y = .value, color = .operator, group = .operator),
        linewidth = po$ewma$line_width, linetype = po$ewma$line_type, alpha = po$ewma$alpha
      )
    } else {
      p <- p + ggplot2::geom_line(
        data = ew, ggplot2::aes(x = .x, y = .value, group = .operator),
        color = po$ewma$color, linewidth = po$ewma$line_width,
        linetype = po$ewma$line_type, alpha = po$ewma$alpha
      )
    }
  }

  # Target line and label.
  if (isTRUE(po$target$show) && !is.null(s$target) && is.finite(s$target)) {
    p <- p + ggplot2::geom_hline(
      yintercept = s$target, color = po$target$color,
      linewidth = po$target$line_width, linetype = po$target$line_type, alpha = po$target$alpha
    )
    if (isTRUE(po$target$label)) {
      tx <- do.call(rbind, lapply(split(raw, raw$.operator), function(z) {
        data.frame(.operator = z$.operator[1L], .x = max(z$.x, na.rm = TRUE), .value = s$target)
      }))
      lab <- .tablearn_or(po$target$label_text, paste0("Target: ", signif(s$target, 4)))
      tx$.label <- lab
      p <- p + ggplot2::geom_text(
        data = tx, ggplot2::aes(x = .x, y = .value, label = .label),
        inherit.aes = FALSE, hjust = 1.02, vjust = -0.45,
        size = po$target$label_size, color = po$target$color
      )
    }
  }

  # Breakpoint lines and labels on the case-number axis.
  if (isTRUE(po$breakpoint$show) && !use_order) {
    bpdf <- do.call(rbind, lapply(names(x$breakpoint), function(nm) {
      bp <- x$breakpoint[[nm]]
      if (is.null(bp) || !is.finite(bp$breakpoint)) return(NULL)
      data.frame(.operator = nm, .bp = bp$breakpoint, stringsAsFactors = FALSE)
    }))
    if (!is.null(bpdf) && nrow(bpdf)) {
      p <- p + ggplot2::geom_vline(
        data = bpdf, ggplot2::aes(xintercept = .bp),
        inherit.aes = FALSE, color = po$breakpoint$color,
        linewidth = po$breakpoint$line_width, linetype = po$breakpoint$line_type,
        alpha = po$breakpoint$alpha
      )
      if (isTRUE(po$breakpoint$label)) {
        bpdf$.label <- if (!is.null(po$breakpoint$label_text)) po$breakpoint$label_text else paste0("Case ", round(bpdf$.bp, 1))
        p <- p + ggplot2::geom_text(
          data = bpdf, ggplot2::aes(x = .bp, y = Inf, label = .label),
          inherit.aes = FALSE, angle = po$breakpoint$label_angle,
          hjust = po$breakpoint$label_hjust, vjust = po$breakpoint$label_vjust,
          size = po$breakpoint$label_size, color = po$breakpoint$color
        )
      }
    }
  }

  # Estimated proficiency line. This is deliberately separate from the change-point
  # line because a change in slope is not automatically equivalent to competency.
  if (isTRUE(po$proficiency$show) && is.data.frame(prof) && nrow(prof) && !use_order) {
    pp <- prof[prof$achieved %in% TRUE & is.finite(prof$proficiency_case), , drop = FALSE]
    if (nrow(pp)) {
      pp$.operator <- pp$operator
      pp$.prof <- pp$proficiency_case
      p <- p + ggplot2::geom_vline(
        data = pp, ggplot2::aes(xintercept = .prof),
        inherit.aes = FALSE, color = po$proficiency$color,
        linewidth = po$proficiency$line_width, linetype = po$proficiency$line_type,
        alpha = po$proficiency$alpha
      )
      if (isTRUE(po$proficiency$label)) {
        if (!is.null(po$proficiency$label_text)) {
          pp$.label <- po$proficiency$label_text
        } else {
          pp$.label <- paste0(
            .tablearn_lang(s$lang, "Proficiency: case ", "Proficiency: ca "),
            round(pp$.prof, 1)
          )
        }
        p <- p + ggplot2::geom_text(
          data = pp, ggplot2::aes(x = .prof, y = Inf, label = .label),
          inherit.aes = FALSE, angle = po$proficiency$label_angle,
          hjust = po$proficiency$label_hjust, vjust = po$proficiency$label_vjust,
          size = po$proficiency$label_size, color = po$proficiency$color
        )
      }
    }
  }

  # Current point and optional label.
  if (isTRUE(po$current$show) && nrow(cur) && nrow(win)) {
    cc <- merge(
      cur, win[, c(".operator", ".end", ".x")],
      by.x = c("operator", "case"), by.y = c(".operator", ".end"), all.x = TRUE
    )
    if (nrow(cc)) {
      p <- p + ggplot2::geom_point(
        data = cc, ggplot2::aes(x = .x, y = estimate),
        color = po$current$color, fill = po$current$fill,
        shape = po$current$shape, size = po$current$size,
        alpha = po$current$alpha, stroke = po$current$stroke
      )
      if (isTRUE(po$current$label)) {
        if (!is.null(po$current$label_text)) {
          cc$.label <- po$current$label_text
        } else if (s$type == "binary") {
          cc$.label <- paste0("Current: ", format(round(100 * cc$estimate, 1), trim = TRUE), "%")
        } else {
          cc$.label <- paste0("Current: ", signif(cc$estimate, 4))
        }
        p <- p + ggplot2::geom_text(
          data = cc, ggplot2::aes(x = .x, y = estimate, label = .label),
          hjust = po$current$hjust, vjust = po$current$vjust,
          size = po$current$label_size, color = po$current$color
        )
      }
    }
  }

  # Faceting for multiple operators.
  if (multi_op && po$operator$display == "facet") {
    p <- p + ggplot2::facet_wrap(
      stats::as.formula("~ .operator"), ncol = po$facet$ncol,
      nrow = po$facet$nrow, scales = po$facet$scales
    )
  }

  # Optional manual operator palettes for overlay mode.
  if (multi_op && po$operator$display == "overlay") {
    if (!is.null(po$operator$colors)) {
      p <- p + ggplot2::scale_color_manual(values = po$operator$colors, name = po$legend$title)
      p <- p + ggplot2::scale_fill_manual(values = po$operator$colors, name = po$legend$title)
    }
    if (!is.null(po$operator$shapes)) p <- p + ggplot2::scale_shape_manual(values = po$operator$shapes, name = po$legend$title)
    if (!is.null(po$operator$line_types)) p <- p + ggplot2::scale_linetype_manual(values = po$operator$line_types, name = po$legend$title)
  }

  yl <- .tablearn_or(ylab, s$outcome)
  xl <- .tablearn_or(
    xlab,
    if (use_order) .tablearn_or(s$order, "Order") else .tablearn_lang(s$lang, "Consecutive case number", "S\u1ed1 ca li\u00ean ti\u1ebfp")
  )
  tt <- .tablearn_or(
    title,
    .tablearn_lang(s$lang, paste0("Learning curve: ", s$outcome), paste0("\u0110\u01b0\u1eddng cong h\u1ecdc t\u1eadp: ", s$outcome))
  )
  p <- p + ggplot2::labs(
    title = tt, subtitle = subtitle, caption = caption,
    x = xl, y = yl, color = po$legend$title, fill = po$legend$title,
    shape = po$legend$title, linetype = po$legend$title
  )

  ax <- po$axes
  ytrans <- if (isTRUE(ax$y_reverse)) "reverse" else .tablearn_or(ax$y_trans, "identity")
  yexpand <- if (is.null(ax$y_expand)) ggplot2::waiver() else ggplot2::expansion(mult = ax$y_expand)
  is_pct <- identical(ax$y_percent, TRUE) || (identical(ax$y_percent, "auto") && s$type == "binary")
  ylabels <- if (is_pct) {
    acc <- suppressWarnings(as.numeric(ax$percent_accuracy)[1L])
    if (!is.finite(acc) || acc <= 0) acc <- 1
    digs <- max(0L, as.integer(ceiling(-log10(acc))))
    function(z) paste0(formatC(100 * z, format = "f", digits = digs), "%")
  } else ggplot2::waiver()
  p <- p + ggplot2::scale_y_continuous(
    labels = ylabels, limits = ax$y_limits, breaks = ax$y_breaks,
    trans = ytrans, expand = yexpand
  )

  if (!use_order) {
    xtrans <- if (isTRUE(ax$x_reverse)) "reverse" else .tablearn_or(ax$x_trans, "identity")
    xexpand <- if (is.null(ax$x_expand)) ggplot2::waiver() else ggplot2::expansion(mult = ax$x_expand)
    p <- p + ggplot2::scale_x_continuous(
      limits = ax$x_limits, breaks = ax$x_breaks,
      trans = xtrans, expand = xexpand
    )
  } else if (inherits(raw$.x, "Date")) {
    xexpand <- if (is.null(ax$x_expand)) ggplot2::waiver() else ggplot2::expansion(mult = ax$x_expand)
    p <- p + ggplot2::scale_x_date(
      limits = ax$x_limits,
      date_labels = .tablearn_or(ax$x_date_format, "%b %Y"),
      date_breaks = .tablearn_or(ax$x_date_breaks, "1 month"),
      expand = xexpand
    )
  }

  p <- p + ggplot2::coord_cartesian(clip = .tablearn_or(ax$clip, "on"))
  p + .tablearn_theme(po)
}

.tablearn_plot_cusum <- function(x, po) {
  pieces <- list()
  signals <- list()

  if (nrow(x$cusum)) {
    d <- x$cusum[, c(".case", ".value", ".operator"), drop = FALSE]
    d$.series <- "CUSUM"
    pieces[[length(pieces) + 1L]] <- d
    if (".signal" %in% names(x$cusum) && any(x$cusum$.signal, na.rm = TRUE)) {
      q <- x$cusum[x$cusum$.signal, c(".case", ".value", ".operator"), drop = FALSE]
      q$.series <- "CUSUM"; signals[[length(signals) + 1L]] <- q
    }
  }
  if (nrow(x$lccusum)) {
    d <- x$lccusum[, c(".case", ".value", ".operator"), drop = FALSE]
    d$.series <- "LC-CUSUM"
    pieces[[length(pieces) + 1L]] <- d
    if (".proficient" %in% names(x$lccusum) && any(x$lccusum$.proficient, na.rm = TRUE)) {
      q <- x$lccusum[x$lccusum$.proficient, c(".case", ".value", ".operator"), drop = FALSE]
      q$.series <- "LC-CUSUM"; signals[[length(signals) + 1L]] <- q
    }
  }
  if (nrow(x$racusum)) {
    d <- x$racusum[, c(".case", ".value", ".operator"), drop = FALSE]
    d$.series <- "RA-CUSUM"
    pieces[[length(pieces) + 1L]] <- d
    if (".signal" %in% names(x$racusum) && any(x$racusum$.signal, na.rm = TRUE)) {
      q <- x$racusum[x$racusum$.signal, c(".case", ".value", ".operator"), drop = FALSE]
      q$.series <- "RA-CUSUM"; signals[[length(signals) + 1L]] <- q
    }
  }
  if (!length(pieces)) return(NULL)

  d <- do.call(rbind, pieces)
  cp <- po$cusum
  p <- ggplot2::ggplot(
    d,
    ggplot2::aes(x = .case, y = .value, color = .series, linetype = .series,
                 group = interaction(.operator, .series))
  )

  if (isTRUE(cp$zero_show)) {
    p <- p + ggplot2::geom_hline(
      yintercept = 0, linewidth = cp$zero_width,
      color = cp$zero_color, linetype = cp$zero_type
    )
  }

  p <- p + ggplot2::geom_line(linewidth = cp$line_width, alpha = cp$alpha)
  if (isTRUE(cp$points)) {
    p <- p + ggplot2::geom_point(size = cp$point_size, shape = cp$point_shape, alpha = cp$alpha)
  }

  # Decision limits: LC-CUSUM may be operator-specific; standard/RA limits are
  # supplied as scalar settings and are drawn globally.
  if (isTRUE(cp$decision_show)) {
    if (nrow(x$lccusum) && ".limit" %in% names(x$lccusum)) {
      ld <- unique(x$lccusum[, c(".operator", ".limit"), drop = FALSE])
      ld <- ld[is.finite(ld$.limit), , drop = FALSE]
      if (nrow(ld)) {
        p <- p + ggplot2::geom_hline(
          data = ld, ggplot2::aes(yintercept = .limit), inherit.aes = FALSE,
          color = cp$decision_color, linewidth = cp$decision_width,
          linetype = cp$decision_type
        )
      }
    }
    if (!is.null(x$settings$cusum_limit) && is.finite(x$settings$cusum_limit)) {
      p <- p + ggplot2::geom_hline(
        yintercept = x$settings$cusum_limit, color = cp$decision_color,
        linewidth = cp$decision_width, linetype = cp$decision_type
      )
    }
    if (!is.null(x$settings$racusum_limit) && is.finite(x$settings$racusum_limit)) {
      p <- p + ggplot2::geom_hline(
        yintercept = x$settings$racusum_limit, color = cp$decision_color,
        linewidth = cp$decision_width, linetype = cp$decision_type
      )
    }
  }

  if (isTRUE(cp$signal_show) && length(signals)) {
    sg <- do.call(rbind, signals)
    p <- p + ggplot2::geom_point(
      data = sg, ggplot2::aes(x = .case, y = .value), inherit.aes = FALSE,
      color = cp$signal_color, fill = cp$signal_fill,
      shape = cp$signal_shape, size = cp$signal_size, stroke = cp$signal_stroke
    )
  }

  if (length(unique(d$.operator)) > 1L && po$operator$display == "facet") {
    p <- p + ggplot2::facet_wrap(
      stats::as.formula("~ .operator"), ncol = po$facet$ncol,
      nrow = po$facet$nrow, scales = po$facet$scales
    )
  }

  series_present <- unique(d$.series)
  cols <- cp$colors
  if (!is.null(cols)) {
    cols <- cols[names(cols) %in% series_present]
    if (length(cols)) p <- p + ggplot2::scale_color_manual(values = cols)
  }
  lts <- cp$line_types
  if (!is.null(lts)) {
    lts <- lts[names(lts) %in% series_present]
    if (length(lts)) p <- p + ggplot2::scale_linetype_manual(values = lts)
  }

  ttl <- .tablearn_or(
    cp$title,
    .tablearn_lang(x$settings$lang, "Sequential performance monitoring", "Gi\u00e1m s\u00e1t hi\u1ec7u su\u1ea5t theo tr\u00ecnh t\u1ef1")
  )
  xlb <- .tablearn_or(cp$xlab, .tablearn_lang(x$settings$lang, "Consecutive case number", "S\u1ed1 ca li\u00ean ti\u1ebfp"))
  ylb <- .tablearn_or(cp$ylab, "CUSUM")
  p <- p + ggplot2::labs(
    title = ttl, subtitle = cp$subtitle, caption = cp$caption,
    x = xlb, y = ylb, color = NULL, linetype = NULL
  )

  p <- p + ggplot2::scale_x_continuous(limits = cp$x_limits, breaks = cp$x_breaks)
  p <- p + ggplot2::scale_y_continuous(limits = cp$y_limits, breaks = cp$y_breaks)
  p + .tablearn_theme(po)
}


# -----------------------------------------------------------------------------
# Publication-ready tables and dependency-light Viewer/base plotting
# -----------------------------------------------------------------------------

.tablearn_fmt_num <- function(x, digits = 2L) {
  x <- suppressWarnings(as.numeric(x))
  out <- rep("", length(x))
  ok <- is.finite(x)
  out[ok] <- formatC(x[ok], format = "f", digits = digits)
  out
}

.tablearn_fmt_p <- function(x, digits = 3L) {
  x <- suppressWarnings(as.numeric(x))
  out <- rep("", length(x))
  ok <- is.finite(x)
  if (!any(ok)) return(out)
  cut <- 10^(-digits)
  out[ok & x < cut] <- paste0("<", formatC(cut, format = "f", digits = digits))
  ii <- ok & x >= cut
  out[ii] <- formatC(x[ii], format = "f", digits = digits)
  out
}

.tablearn_fmt_ci <- function(est, lower, upper, digits = 2L, percent = FALSE) {
  est <- suppressWarnings(as.numeric(est))
  lower <- suppressWarnings(as.numeric(lower))
  upper <- suppressWarnings(as.numeric(upper))
  if (isTRUE(percent)) {
    est <- 100 * est
    lower <- 100 * lower
    upper <- 100 * upper
  }
  e <- .tablearn_fmt_num(est, digits)
  lo <- .tablearn_fmt_num(lower, digits)
  hi <- .tablearn_fmt_num(upper, digits)
  if (isTRUE(percent)) {
    e[nzchar(e)] <- paste0(e[nzchar(e)], "%")
    lo[nzchar(lo)] <- paste0(lo[nzchar(lo)], "%")
    hi[nzchar(hi)] <- paste0(hi[nzchar(hi)], "%")
  }
  out <- e
  have_ci <- nzchar(lo) & nzchar(hi)
  out[have_ci] <- paste0(e[have_ci], " (", lo[have_ci], " to ", hi[have_ci], ")")
  out
}

.tablearn_yes_no <- function(x) {
  out <- rep("", length(x))
  out[!is.na(x) & x] <- "Yes"
  out[!is.na(x) & !x] <- "No"
  out
}

.tablearn_breakpoint_table <- function(x, digit = 2L, p_digit = 3L) {
  if (!length(x$breakpoint)) return(data.frame())
  rows <- lapply(names(x$breakpoint), function(nm) {
    z <- x$breakpoint[[nm]]
    if (is.null(z)) return(NULL)
    ci <- z$ci
    ci_txt <- if (length(ci) >= 2L && all(is.finite(ci[1:2]))) {
      paste0(.tablearn_fmt_num(ci[1], digit), " to ", .tablearn_fmt_num(ci[2], digit))
    } else ""
    data.frame(
      Operator = nm,
      `Change point` = .tablearn_fmt_num(z$breakpoint, digit),
      `95% CI` = ci_txt,
      Detected = if (isTRUE(z$detected)) "Yes" else "No",
      `Slope before` = .tablearn_fmt_num(z$slope_before, digit),
      `Slope after` = .tablearn_fmt_num(z$slope_after, digit),
      `Delta BIC` = .tablearn_fmt_num(z$delta_bic, digit),
      p = .tablearn_fmt_p(z$p_change, p_digit),
      stringsAsFactors = FALSE, check.names = FALSE
    )
  })
  rows <- rows[!vapply(rows, is.null, logical(1))]
  if (!length(rows)) return(data.frame())
  do.call(rbind, rows)
}

.tablearn_monitoring_summary <- function(x, digit = 2L) {
  one <- function(dat, label, limit_col = NULL, signal_col = NULL, forced_limit = NULL) {
    if (!is.data.frame(dat) || !nrow(dat)) return(NULL)
    spl <- split(dat, dat$.operator)
    do.call(rbind, lapply(names(spl), function(op) {
      z <- spl[[op]]
      sig <- NA_real_
      if (!is.null(signal_col) && signal_col %in% names(z)) {
        ii <- which(!is.na(z[[signal_col]]) & as.logical(z[[signal_col]]))
        if (length(ii)) sig <- z$.case[ii[1L]]
      }
      lim <- if (!is.null(forced_limit) && length(forced_limit) == 1L && is.finite(forced_limit)) {
        as.numeric(forced_limit)
      } else if (!is.null(limit_col) && limit_col %in% names(z)) {
        v <- z[[limit_col]][is.finite(z[[limit_col]])]
        if (length(v)) v[1L] else NA_real_
      } else NA_real_
      data.frame(
        Operator = op,
        Method = label,
        `Last case` = max(z$.case, na.rm = TRUE),
        `Current value` = .tablearn_fmt_num(tail(z$.value, 1L), digit),
        `Decision limit` = .tablearn_fmt_num(lim, digit),
        `First signal case` = .tablearn_fmt_num(sig, 0L),
        stringsAsFactors = FALSE, check.names = FALSE
      )
    }))
  }
  rows <- list(
    one(x$cusum, "CUSUM", signal_col = ".signal", forced_limit = x$settings$cusum_limit),
    one(x$lccusum, "LC-CUSUM", limit_col = ".limit", signal_col = ".proficient"),
    one(x$racusum, "RA-CUSUM", signal_col = ".signal", forced_limit = x$settings$racusum_limit)
  )
  rows <- rows[!vapply(rows, is.null, logical(1))]
  if (!length(rows)) return(data.frame())
  out <- do.call(rbind, rows)
  rownames(out) <- NULL
  out
}

.tablearn_make_tables <- function(x, digit = 2L, p_digit = 3L) {
  is_binary <- identical(x$settings$type, "binary")
  nops <- length(unique(x$raw$.operator))
  overview <- data.frame(
    Statistic = c(
      "Outcome", "Outcome type", "Analyzed observations", "Order",
      "Operator", "Window", "Window summary", "Better performance",
      "Target", "Change-point inference"
    ),
    Value = c(
      x$settings$outcome,
      x$settings$type,
      as.character(nrow(x$raw)),
      if (is.null(x$settings$order)) "Current row order" else x$settings$order,
      if (is.null(x$settings$operator)) "Overall" else paste0(x$settings$operator, " (", nops, " operators)"),
      paste0(x$settings$window_type, ", n = ", x$settings$window,
             ", step = ", x$settings$window_step, ", align = ", x$settings$window_align),
      if (identical(x$settings$window_fun, "auto")) if (is_binary) "proportion" else "mean" else x$settings$window_fun,
      x$settings$better,
      if (is.null(x$settings$target)) {
        "Not specified"
      } else if (is_binary) {
        paste0(.tablearn_fmt_num(100 * x$settings$target, digit), "%")
      } else {
        .tablearn_fmt_num(x$settings$target, digit)
      },
      if (isTRUE(x$settings$breakpoint)) paste0("Piecewise model on ", x$settings$break_on, " cases") else "Not requested"
    ),
    stringsAsFactors = FALSE, check.names = FALSE
  )

  current <- data.frame()
  if (is.data.frame(x$current) && nrow(x$current)) {
    z <- x$current
    current <- data.frame(
      Operator = z$operator,
      Case = z$case,
      Window = paste0(z$window_start, "-", z$window_end),
      N = z$n,
      stringsAsFactors = FALSE, check.names = FALSE
    )
    if (is_binary && "events" %in% names(z)) current$Events <- z$events
    current$`Estimate (95% CI)` <- .tablearn_fmt_ci(z$estimate, z$lower, z$upper, digit, percent = is_binary)
    if (is_binary) {
      current$`Change (percentage points)` <- .tablearn_fmt_num(100 * z$change, digit)
    } else {
      current$Change <- .tablearn_fmt_num(z$change, digit)
    }
    if (any(is.finite(z$target))) {
      current$Target <- if (is_binary) {
        paste0(.tablearn_fmt_num(100 * z$target, digit), "%")
      } else {
        .tablearn_fmt_num(z$target, digit)
      }
      current$`Target met` <- .tablearn_yes_no(z$target_met)
    }
  }

  proficiency <- data.frame()
  if (is.data.frame(x$proficiency) && nrow(x$proficiency)) {
    z <- x$proficiency
    proficiency <- data.frame(
      Operator = z$operator,
      Status = z$status,
      `Proficiency case` = .tablearn_fmt_num(z$proficiency_case, 0L),
      Method = z$method,
      Evidence = z$evidence,
      `Current phase` = z$current_phase,
      `Cases since proficiency` = .tablearn_fmt_num(z$cases_since_proficiency, 0L),
      stringsAsFactors = FALSE, check.names = FALSE
    )
  }

  evidence <- data.frame()
  if (is.data.frame(x$proficiency_evidence) && nrow(x$proficiency_evidence)) {
    z <- x$proficiency_evidence
    evidence <- data.frame(
      Operator = z$operator,
      Criterion = z$criterion,
      Available = .tablearn_yes_no(z$available),
      Met = .tablearn_yes_no(z$met),
      Case = .tablearn_fmt_num(z$case, 0L),
      Detail = z$detail,
      stringsAsFactors = FALSE, check.names = FALSE
    )
  }

  phases <- data.frame()
  if (is.data.frame(x$phases) && nrow(x$phases)) {
    z <- x$phases
    phases <- data.frame(
      Operator = z$operator,
      Phase = z$phase,
      Cases = paste0(z$start, "-", z$end),
      N = z$n,
      stringsAsFactors = FALSE, check.names = FALSE
    )
    if (is_binary && "events" %in% names(z)) phases$Events <- z$events
    phases$`Estimate (95% CI)` <- .tablearn_fmt_ci(z$estimate, z$lower, z$upper, digit, percent = is_binary)
  }

  windows <- data.frame()
  if (is.data.frame(x$window) && nrow(x$window)) {
    z <- x$window
    windows <- data.frame(
      Operator = z$.operator,
      Window = z$.window_label,
      N = z$.n,
      stringsAsFactors = FALSE, check.names = FALSE
    )
    if (is_binary && ".events" %in% names(z)) windows$Events <- z$.events
    windows$`Estimate (95% CI)` <- .tablearn_fmt_ci(z$.value, z$.lower, z$.upper, digit, percent = is_binary)
  }

  sensitivity <- data.frame()
  if (is.data.frame(x$sensitivity) && nrow(x$sensitivity)) {
    z <- x$sensitivity
    sensitivity <- data.frame(
      Operator = z$operator,
      `Window n` = z$window,
      Points = z$points,
      `Current estimate` = if (is_binary) paste0(.tablearn_fmt_num(100 * z$current, digit), "%") else .tablearn_fmt_num(z$current, digit),
      `Change point` = .tablearn_fmt_num(z$breakpoint, digit),
      Detected = .tablearn_yes_no(z$detected),
      stringsAsFactors = FALSE, check.names = FALSE
    )
  }

  list(
    Overview = overview,
    Current_performance = current,
    Change_point = .tablearn_breakpoint_table(x, digit, p_digit),
    Proficiency = proficiency,
    Proficiency_evidence = evidence,
    Phase_classification = phases,
    Window_performance = windows,
    Window_sensitivity = sensitivity,
    Sequential_monitoring = .tablearn_monitoring_summary(x, digit)
  )
}

.tablearn_has_cusum <- function(x) {
  (is.data.frame(x$cusum) && nrow(x$cusum)) ||
    (is.data.frame(x$lccusum) && nrow(x$lccusum)) ||
    (is.data.frame(x$racusum) && nrow(x$racusum))
}

.tablearn_plot_x <- function(raw, cases, use_order = FALSE) {
  cases <- as.numeric(cases)
  if (!isTRUE(use_order)) return(cases)
  rr <- raw[order(raw$.case), , drop = FALSE]
  if (!nrow(rr)) return(cases)
  if (inherits(rr$.order, "Date")) {
    val <- stats::approx(rr$.case, as.numeric(rr$.order), xout = cases, rule = 2)$y
    return(as.Date(val, origin = "1970-01-01"))
  }
  if (is.numeric(rr$.order)) {
    return(stats::approx(rr$.case, rr$.order, xout = cases, rule = 2)$y)
  }
  cases
}

.tablearn_base_learning_panel <- function(x, op = NULL, overlay = FALSE, op_colors = NULL) {
  po <- x$plot_options
  s <- x$settings
  raw <- x$raw
  use_order <- identical(s$x_axis, "order") && (inherits(raw$.order, "Date") || is.numeric(raw$.order))
  if (!is.null(op)) raw <- raw[raw$.operator == op, , drop = FALSE]
  if (!nrow(raw)) return(invisible(NULL))

  get_part <- function(dat) {
    if (!is.data.frame(dat) || !nrow(dat)) return(dat)
    if (is.null(op) || isTRUE(overlay)) dat else dat[dat$.operator == op, , drop = FALSE]
  }
  win <- get_part(x$window); sm <- get_part(x$smooth); ew <- get_part(x$ewma)
  cur <- get_part(x$current); ph <- get_part(x$phases); prof <- get_part(x$proficiency)

  xraw <- if (use_order) raw$.order else raw$.case
  yvals <- raw$.value
  for (z in list(win$.value, win$.lower, win$.upper, sm$.value, sm$.lower, sm$.upper, ew$.value)) {
    if (length(z)) yvals <- c(yvals, z)
  }
  if (!is.null(s$target) && is.finite(s$target)) yvals <- c(yvals, s$target)
  yvals <- yvals[is.finite(yvals)]
  if (!length(yvals)) yvals <- c(0, 1)
  yr <- range(yvals)
  if (identical(s$type, "binary")) yr <- c(0, 1)
  if (diff(yr) <= 0) yr <- yr + c(-0.5, 0.5)
  pad <- diff(yr) * 0.06
  yr <- yr + c(-pad, pad)
  if (identical(s$type, "binary")) yr <- c(0, 1)

  xlab <- if (!is.null(s$order) && use_order) s$order else "Consecutive case number"
  ylab <- if (identical(s$type, "binary")) "Event rate (%)" else s$outcome
  ttl <- .tablearn_or(s$title, paste0("Learning curve: ", s$outcome))
  if (!is.null(op) && !identical(op, "Overall")) ttl <- paste0(ttl, " \u2014 ", op)

  graphics::plot(xraw, raw$.value, type = "n", xlab = xlab, ylab = ylab,
                 ylim = yr, main = ttl, las = 1,
                 yaxt = if (identical(s$type, "binary")) "n" else "s")
  if (identical(s$type, "binary")) {
    at <- seq(0, 1, by = 0.25)
    graphics::axis(2, at = at, labels = paste0(round(100 * at), "%"), las = 1)
  }
  graphics::grid(col = "#E5E7EB", lty = 1)

  # Phase shading is intentionally used only for a single/faceted operator.
  if (isTRUE(po$phase$show) && nrow(ph) && !isTRUE(overlay)) {
    fills <- .tablearn_or(po$phase$fills, c("#F5B7B1", "#F9E79F", "#ABEBC6"))
    for (i in seq_len(nrow(ph))) {
      xx <- .tablearn_plot_x(raw, c(ph$start[i] - 0.5, ph$end[i] + 0.5), use_order)
      col <- grDevices::adjustcolor(fills[(i - 1L) %% length(fills) + 1L], alpha.f = .tablearn_or(po$phase$alpha, .10))
      graphics::rect(xx[1], graphics::par("usr")[3], xx[2], graphics::par("usr")[4], col = col, border = NA)
    }
    graphics::box()
  }

  ops <- unique(raw$.operator)
  if (isTRUE(overlay) && length(ops) > 1L) {
    if (is.null(op_colors) || length(op_colors) < length(ops)) op_colors <- grDevices::hcl.colors(length(ops), "Dark 3")
    for (i in seq_along(ops)) {
      oo <- ops[i]
      rr <- raw[raw$.operator == oo, , drop = FALSE]
      ww <- if (nrow(win)) win[win$.operator == oo, , drop = FALSE] else win
      xr <- if (use_order) rr$.order else rr$.case
      graphics::points(xr, rr$.value, pch = 16, cex = .65,
                       col = grDevices::adjustcolor(op_colors[i], alpha.f = .28))
      if (nrow(ww)) {
        xw <- if (use_order) ww$.x_order else ww$.x_case
        graphics::lines(xw, ww$.value, col = op_colors[i], lwd = 2)
      }
    }
    graphics::legend("topright", legend = ops, col = op_colors, lwd = 2, bty = "n", cex = .85)
  } else {
    if (isTRUE(po$raw$show)) {
      graphics::points(xraw, raw$.value, pch = 16, cex = .65,
                       col = grDevices::adjustcolor(.tablearn_or(po$raw$color, "#7A7A7A"),
                                                   alpha.f = .tablearn_or(po$raw$alpha, .28)))
    }
    if (nrow(win)) {
      xw <- if (use_order) win$.x_order else win$.x_case
      if (isTRUE(po$window_interval$show) && any(is.finite(win$.lower) & is.finite(win$.upper))) {
        ok <- is.finite(win$.lower) & is.finite(win$.upper)
        graphics::segments(xw[ok], win$.lower[ok], xw[ok], win$.upper[ok],
                           col = grDevices::adjustcolor(.tablearn_or(po$window_interval$color, "#1F5A94"), alpha.f = .45))
      }
      if (isTRUE(po$window$show)) {
        graphics::lines(xw, win$.value, col = .tablearn_or(po$window$line_color, "#1F5A94"),
                        lwd = max(1, .tablearn_or(po$window$line_width, 1.05) * 1.5),
                        lty = .tablearn_or(po$window$line_type, 1))
        graphics::points(xw, win$.value, pch = 21, bg = .tablearn_or(po$window$fill, "white"),
                         col = .tablearn_or(po$window$color, "#1F5A94"), cex = .78)
      }
    }
    if (nrow(sm) && isTRUE(po$smooth$show)) {
      xs <- .tablearn_plot_x(raw, sm$.case, use_order)
      if (isTRUE(po$smooth$ci_show) && any(is.finite(sm$.lower) & is.finite(sm$.upper))) {
        ok <- is.finite(sm$.lower) & is.finite(sm$.upper)
        if (sum(ok) >= 2L) {
          graphics::polygon(c(xs[ok], rev(xs[ok])), c(sm$.lower[ok], rev(sm$.upper[ok])),
                            col = grDevices::adjustcolor(.tablearn_or(po$smooth$fill, "#B23A48"),
                                                        alpha.f = .tablearn_or(po$smooth$ci_alpha, .10)),
                            border = NA)
        }
      }
      graphics::lines(xs, sm$.value, col = .tablearn_or(po$smooth$color, "#B23A48"),
                      lwd = max(1, .tablearn_or(po$smooth$line_width, 1.2) * 1.5),
                      lty = .tablearn_or(po$smooth$line_type, 1))
    }
    if (nrow(ew) && isTRUE(po$ewma$show)) {
      xe <- if (use_order) ew$.order else ew$.case
      graphics::lines(xe, ew$.value, col = .tablearn_or(po$ewma$color, "#7A3E9D"),
                      lwd = 1.5, lty = .tablearn_or(po$ewma$line_type, 2))
    }
  }

  if (!is.null(s$target) && is.finite(s$target) && isTRUE(po$target$show)) {
    graphics::abline(h = s$target, col = .tablearn_or(po$target$color, "#2E7D32"),
                     lty = .tablearn_or(po$target$line_type, 3), lwd = 1.4)
  }
  if (!isTRUE(overlay)) {
    bp <- x$breakpoint[[if (is.null(op)) names(x$breakpoint)[1L] else op]]
    if (!is.null(bp) && isTRUE(bp$detected) && isTRUE(po$breakpoint$show)) {
      xb <- .tablearn_plot_x(raw, bp$breakpoint, use_order)
      graphics::abline(v = xb, col = .tablearn_or(po$breakpoint$color, "#222222"),
                       lty = .tablearn_or(po$breakpoint$line_type, 2), lwd = 1.2)
    }
    if (nrow(prof) && isTRUE(prof$achieved[1L]) && is.finite(prof$proficiency_case[1L]) && isTRUE(po$proficiency$show)) {
      xp <- .tablearn_plot_x(raw, prof$proficiency_case[1L], use_order)
      graphics::abline(v = xp, col = .tablearn_or(po$proficiency$color, "#00796B"),
                       lty = .tablearn_or(po$proficiency$line_type, 1), lwd = 1.5)
    }
    if (nrow(cur) && isTRUE(po$current$show)) {
      xc <- .tablearn_plot_x(raw, cur$case[1L], use_order)
      graphics::points(xc, cur$estimate[1L], pch = 21, bg = .tablearn_or(po$current$fill, "#FFD166"),
                       col = .tablearn_or(po$current$color, "#000000"), cex = 1.25, lwd = 1)
    }
  }
  invisible(NULL)
}

.tablearn_plot_base_learning <- function(x) {
  ops <- unique(x$raw$.operator)
  overlay <- length(ops) > 1L && identical(x$settings$operator_display, "overlay")
  if (overlay) return(.tablearn_base_learning_panel(x, overlay = TRUE))

  old <- graphics::par(no.readonly = TRUE)
  on.exit(graphics::par(old), add = TRUE)
  if (length(ops) > 1L) {
    nc <- ceiling(sqrt(length(ops)))
    nr <- ceiling(length(ops) / nc)
    graphics::par(mfrow = c(nr, nc), mar = c(4.2, 4.2, 3.4, 1.1))
  }
  for (op in ops) .tablearn_base_learning_panel(x, op = op, overlay = FALSE)
  invisible(x)
}

.tablearn_plot_base_cusum <- function(x) {
  if (!.tablearn_has_cusum(x)) .tablearn_stop("No CUSUM analysis is available in this result.")
  cp <- x$plot_options$cusum
  dats <- list(CUSUM = x$cusum, `LC-CUSUM` = x$lccusum, `RA-CUSUM` = x$racusum)
  dats <- dats[vapply(dats, function(z) is.data.frame(z) && nrow(z), logical(1))]
  ops <- unique(unlist(lapply(dats, function(z) z$.operator), use.names = FALSE))
  old <- graphics::par(no.readonly = TRUE)
  on.exit(graphics::par(old), add = TRUE)
  if (length(ops) > 1L) {
    nc <- ceiling(sqrt(length(ops))); nr <- ceiling(length(ops) / nc)
    graphics::par(mfrow = c(nr, nc), mar = c(4.2, 4.2, 3.4, 1.1))
  }
  cols <- .tablearn_or(cp$colors, c(CUSUM = "#1F5A94", `LC-CUSUM` = "#B23A48", `RA-CUSUM` = "#7A3E9D"))
  for (op in ops) {
    parts <- lapply(dats, function(z) z[z$.operator == op, , drop = FALSE])
    yr <- range(unlist(lapply(parts, function(z) z$.value)), finite = TRUE)
    lims <- unlist(lapply(parts, function(z) if (".limit" %in% names(z)) z$.limit else numeric()), use.names = FALSE)
    lims <- c(lims, x$settings$cusum_limit, x$settings$racusum_limit)
    lims <- suppressWarnings(as.numeric(lims))
    lims <- lims[is.finite(lims)]
    if (length(lims)) yr <- range(c(yr, lims), finite = TRUE)
    yr <- range(c(yr, 0), finite = TRUE)
    if (!all(is.finite(yr)) || diff(yr) == 0) yr <- c(-1, 1)
    first <- parts[[which(vapply(parts, nrow, integer(1)) > 0L)[1L]]]
    ttl <- .tablearn_or(cp$title, "Sequential performance monitoring")
    if (length(ops) > 1L) ttl <- paste0(ttl, " \u2014 ", op)
    graphics::plot(first$.case, first$.value, type = "n", xlab = .tablearn_or(cp$xlab, "Consecutive case number"),
                   ylab = .tablearn_or(cp$ylab, "CUSUM"), ylim = yr, main = ttl, las = 1)
    graphics::grid(col = "#E5E7EB", lty = 1)
    if (isTRUE(cp$zero_show)) graphics::abline(h = 0, col = .tablearn_or(cp$zero_color, "#777777"), lty = .tablearn_or(cp$zero_type, 1))
    shown <- character()
    for (nm in names(parts)) {
      z <- parts[[nm]]
      if (!nrow(z)) next
      cc <- if (!is.null(cols[[nm]])) cols[[nm]] else "#1F5A94"
      graphics::lines(z$.case, z$.value, col = cc, lwd = max(1.2, .tablearn_or(cp$line_width, 1.05) * 1.5))
      shown <- c(shown, nm)
      if (isTRUE(cp$decision_show)) {
        lim <- if (".limit" %in% names(z)) z$.limit[is.finite(z$.limit)] else numeric()
        if (!length(lim) && identical(nm, "CUSUM") && !is.null(x$settings$cusum_limit) && is.finite(x$settings$cusum_limit)) {
          lim <- x$settings$cusum_limit
        }
        if (!length(lim) && identical(nm, "RA-CUSUM") && !is.null(x$settings$racusum_limit) && is.finite(x$settings$racusum_limit)) {
          lim <- x$settings$racusum_limit
        }
        if (length(lim)) graphics::abline(h = lim[1L], col = .tablearn_or(cp$decision_color, "#222222"), lty = .tablearn_or(cp$decision_type, 2))
      }
      sigcol <- if (nm == "LC-CUSUM") ".proficient" else ".signal"
      if (isTRUE(cp$signal_show) && sigcol %in% names(z)) {
        ii <- which(!is.na(z[[sigcol]]) & z[[sigcol]])
        if (length(ii)) graphics::points(z$.case[ii[1L]], z$.value[ii[1L]], pch = 21,
                                         bg = .tablearn_or(cp$signal_fill, "#FFD166"), col = .tablearn_or(cp$signal_color, "#000000"), cex = 1.2)
      }
    }
    if (length(shown) > 1L) graphics::legend("topright", shown, col = unname(cols[shown]), lwd = 2, bty = "n", cex = .85)
  }
  invisible(x)
}

.tablearn_base64 <- function(bytes) {
  bytes <- as.integer(bytes)
  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 = "")
}

.tablearn_view_plot_html <- function(x, type = c("learning", "cusum"), format = c("png", "svg"),
                                     width = 1200L, height = 780L, res = 144L) {
  type <- match.arg(type); format <- match.arg(format)
  path <- tempfile(fileext = paste0(".", format))
  on.exit(unlink(path), add = TRUE)
  if (format == "png") {
    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)
  } else {
    grDevices::svg(path, width = width / res, height = height / res, onefile = TRUE, bg = "white", family = "sans")
  }
  ok <- TRUE
  tryCatch({
    saved_plot <- x$plots[[type]]
    if (!is.null(saved_plot)) {
      print(saved_plot)
    } else if (type == "learning") {
      .tablearn_plot_base_learning(x)
    } else {
      .tablearn_plot_base_cusum(x)
    }
  }, error = function(e) {
    ok <<- FALSE
  }, finally = {
    try(grDevices::dev.off(), silent = TRUE)
  })
  if (!ok || !file.exists(path) || !is.finite(file.info(path)$size) || file.info(path)$size <= 0) return("")
  raw <- readBin(path, what = "raw", n = file.info(path)$size)
  mime <- if (format == "png") "image/png" else "image/svg+xml"
  alt <- if (type == "learning") "Learning curve" else "CUSUM monitoring"
  paste0('<div class="tablearn-figure"><img alt="', .r4vn_view_escape(alt),
         '" src="data:', mime, ';base64,', .tablearn_base64(raw), '"></div>')
}

.tablearn_view_one <- function(x) {
  t <- x$tables
  blocks <- character()
  if (is.data.frame(t$Overview) && nrow(t$Overview)) {
    blocks <- c(blocks, paste0('<section class="r4vn-section"><h2>Overview</h2>',
                               .r4vn_view_key_values(t$Overview), '</section>'))
  }
  add <- function(title, obj, description = NULL) {
    if (is.data.frame(obj) && nrow(obj)) blocks <<- c(blocks, .r4vn_view_section(title, obj, description = description))
  }
  add("Current performance", t$Current_performance,
      "The latest displayed window summarizes the most recent performance under the selected window definition.")
  add("Change-point analysis", t$Change_point,
      "Piecewise regression is fitted to the selected inferential series; a change point is not automatically equivalent to proficiency.")
  add("Proficiency assessment", t$Proficiency)
  add("Evidence used for proficiency", t$Proficiency_evidence)
  add("Phase classification", t$Phase_classification)
  add("Window-size sensitivity", t$Window_sensitivity)
  add("Sequential monitoring summary", t$Sequential_monitoring)
  add("Window performance", t$Window_performance,
      "Rolling windows may overlap; use them for performance display rather than as independent observations for inference.")

  if (isTRUE(x$settings$interpretation) && nzchar(x$interpretation)) {
    txt <- strsplit(x$interpretation, "\n", fixed = TRUE)[[1L]]
    interp <- data.frame(Interpretation = txt[nzchar(txt)], stringsAsFactors = FALSE, check.names = FALSE)
    add("Interpretation", interp)
  }

  if (isTRUE(x$settings$plot)) {
    fmt <- .tablearn_or(x$settings$viewer_plot_format, "png")
    if (x$settings$plot_type %in% c("learning", "both")) {
      img <- .tablearn_view_plot_html(x, "learning", fmt)
      if (nzchar(img)) blocks <- c(blocks, paste0('<section class="r4vn-section"><h2>Learning curve</h2>', img, '</section>'))
    }
    if (x$settings$plot_type %in% c("cusum", "both") && .tablearn_has_cusum(x)) {
      img <- .tablearn_view_plot_html(x, "cusum", fmt)
      if (nzchar(img)) blocks <- c(blocks, paste0('<section class="r4vn-section"><h2>Sequential monitoring</h2>', img, '</section>'))
    }
  }
  paste0(blocks, collapse = "")
}

.r4vn_viewer_tablearn <- function(x) {
  multi <- inherits(x, "r4vn_tablearn_multi")
  if (multi) {
    body <- paste0(vapply(x$outcome_names, function(nm) {
      z <- x$outcomes[[nm]]
      paste0('<div class="tablearn-outcome"><div class="tablearn-outcome-title">Outcome: ',
             .r4vn_view_escape(nm), '</div>', .tablearn_view_one(z), '</div>')
    }, character(1)), collapse = "")
    ttl <- .tablearn_or(x$settings$title, "Learning-curve analysis")
    subtitle <- paste0("Multiple outcomes: ", paste(x$outcome_names, collapse = ", "))
  } else {
    body <- .tablearn_view_one(x)
    ttl <- .tablearn_or(x$settings$title, paste0("Learning-curve analysis: ", x$settings$outcome))
    subtitle <- "Sequential clinical/procedural performance"
  }
  notes <- c(
    "Standard analysis and Viewer plots use base R; no plotting add-on package is required.",
    "Overlapping rolling windows share observations and should not be treated as independent inferential observations.",
    "Proficiency is distinct from a statistical change point and should be interpreted against a clinically meaningful competency criterion whenever possible."
  )
  doc <- .r4vn_view_document(ttl, body, notes = notes, subtitle = subtitle, prefix = "r4vn-tablearn-")
  extra_css <- paste0(
    '<style>.tablearn-figure{text-align:center;overflow:auto}.tablearn-figure img{max-width:100%;height:auto;border:1px solid #eef2f7;border-radius:8px}',
    '.tablearn-outcome{margin-bottom:26px}.tablearn-outcome-title{font-size:19px;font-weight:700;color:#0f172a;margin:4px 0 12px;padding:10px 2px;border-bottom:2px solid #e2e8f0}</style>'
  )
  doc$html <- sub("</head>", paste0(extra_css, "</head>"), doc$html, fixed = TRUE)
  writeLines(doc$html, doc$file, useBytes = TRUE)
  doc
}

#' @export
print.r4vn_tablearn <- function(x, ...) {
  cat("R4VN learning-curve analysis\n")
  cat("Outcome:", x$settings$outcome, " | Type:", x$settings$type, "\n")
  cat("Window:", x$settings$window_type, "n =", x$settings$window,
      "| step =", x$settings$window_step, "| align =", x$settings$window_align, "\n")
  if (is.list(x$tables) && is.data.frame(x$tables$Current_performance) && nrow(x$tables$Current_performance)) {
    cat("\nCurrent performance\n")
    print(x$tables$Current_performance, row.names = FALSE, right = TRUE)
  }
  if (is.list(x$tables) && is.data.frame(x$tables$Change_point) && nrow(x$tables$Change_point)) {
    cat("\nChange-point analysis\n")
    print(x$tables$Change_point, row.names = FALSE, right = TRUE)
  }
  if (is.list(x$tables) && is.data.frame(x$tables$Proficiency) && nrow(x$tables$Proficiency)) {
    cat("\nProficiency assessment\n")
    print(x$tables$Proficiency, row.names = FALSE, right = TRUE)
  }
  if (is.list(x$tables) && is.data.frame(x$tables$Phase_classification) && nrow(x$tables$Phase_classification)) {
    cat("\nPhase classification\n")
    print(x$tables$Phase_classification, row.names = FALSE, right = TRUE)
  }
  if (isTRUE(x$settings$interpretation)) {
    cat("\nInterpretation\n", x$interpretation, "\n", sep = "")
  }
  invisible(x)
}

#' @export
summary.r4vn_tablearn <- function(object, ...) {
  list(
    settings = object$settings,
    current = object$current,
    current_status = object$current_status,
    breakpoint = object$breakpoint,
    proficiency = object$proficiency,
    proficiency_evidence = object$proficiency_evidence,
    phases = object$phases,
    phase_classification = object$phase_classification,
    sensitivity = object$sensitivity,
    tables = object$tables,
    interpretation = object$interpretation
  )
}

#' @export
plot.r4vn_tablearn <- function(x, type = c("learning", "cusum"), ...) {
  type <- match.arg(type)
  if (identical(type, "cusum") && !.tablearn_has_cusum(x)) {
    .tablearn_stop("No CUSUM result is available. Enable `cusum`, `lccusum`, or `racusum` first.")
  }
  p <- x$plots[[type]]
  if (!is.null(p)) {
    print(p)
    return(invisible(p))
  }
  if (identical(type, "learning")) .tablearn_plot_base_learning(x) else .tablearn_plot_base_cusum(x)
  invisible(x)
}

#' @export
print.r4vn_tablearn_multi <- function(x, ...) {
  cat("R4VN learning-curve analysis: multiple outcomes\n")
  cat("Outcomes:", paste(x$outcome_names, collapse = ", "), "\n")
  for (nm in x$outcome_names) {
    cat("\n--- ", nm, " ---\n", sep = "")
    z <- x$outcomes[[nm]]
    if (nrow(z$current)) print(z$current, row.names = FALSE)
    if (is.data.frame(z$proficiency) && nrow(z$proficiency)) {
      print(z$proficiency[, c("operator", "status", "proficiency_case", "evidence"), drop = FALSE], row.names = FALSE)
    }
  }
  invisible(x)
}