R/tabscore.R

Defines functions .r4vn_score_bootstrap_validation .r4vn_score_show_html .r4vn_score_plot_html .r4vn_score_plot_svg .r4vn_score_plot_png .r4vn_score_base64 .r4vn_score_plot_draw .r4vn_score_plot_title .r4vn_score_available_plots .r4vn_score_view_transpose .r4vn_score_html_table .r4vn_score_html_escape .r4vn_score_risk_display .r4vn_score_scorecard_display .r4vn_score_orientation_display .r4vn_score_calibration_data .r4vn_score_decision_curve .r4vn_score_poisson_perf .r4vn_score_comparison_cox .r4vn_score_comparison_binary .r4vn_score_auc_compare .r4vn_score_refprob_table .r4vn_score_reference_probability .r4vn_score_ipcw_brier .r4vn_score_cut_metrics_cox .r4vn_score_cutoff_table_cox .r4vn_score_time_roc .r4vn_score_km_value .r4vn_score_km_censor .r4vn_score_risk_table_cox .r4vn_score_cox_risk .r4vn_score_risk_table_poisson .r4vn_score_risk_table_binary .r4vn_score_add_observed_poisson .r4vn_score_add_observed_binary .r4vn_score_risk_group .r4vn_score_cutoff_table_binary .r4vn_score_find_risk_score .r4vn_score_roc_table .r4vn_score_cut_metrics .r4vn_score_binom_ci .r4vn_score_model_publication .r4vn_score_model_table .r4vn_score_binary_performance .r4vn_score_calibration .r4vn_score_recal_binary .r4vn_score_theoretical_range .r4vn_score_attainable .r4vn_score_build_points .r4vn_score_manual_points .r4vn_score_match_categories .r4vn_score_normalize_category_key .r4vn_score_apply_dictionary .r4vn_score_perf_metric .r4vn_score_cindex .r4vn_score_auc_ci .r4vn_score_auc .r4vn_score_effect_dictionary .r4vn_score_predict_lp .r4vn_score_validate_scoreref .r4vn_score_scoreref_value .r4vn_score_effect_measure .r4vn_score_reference_row .r4vn_score_transform_apply .r4vn_score_transform_fit .r4vn_score_cut_labels .r4vn_score_easy_cuts .r4vn_score_nice_step .r4vn_score_select .r4vn_score_select_lasso .r4vn_score_select_purposeful .r4vn_score_coef_change .r4vn_score_term_p .r4vn_score_lrt .r4vn_score_fit .r4vn_score_formula .r4vn_score_is_binary .r4vn_score_binary .r4vn_score_rebuild_predictor_tokens .r4vn_score_apply_model_input .r4vn_score_capture_model_input .r4vn_score_prepare_declared .r4vn_score_predictor_names .r4vn_score_predictor_spec .r4vn_score_parse_token .r4vn_score_resolve_name .r4vn_score_extract_model .r4vn_score_active_data .r4vn_score_label .r4vn_score_ci_text .r4vn_score_pct .r4vn_score_num .r4vn_score_bt .r4vn_score_null

# tabscore.R -----------------------------------------------------------------
# R4VN: development, simplification, validation and presentation of clinical
# prediction scorecards.
#
# Risk-only scoring revision: by default, protective model contrasts are
# re-referenced to the lowest-risk category so bedside points are add-only
# (all scoring OR/HR/RR/IRR >= 1 and all clinical points >= 0).
#
# This file is intentionally self-contained. Optional packages are used only
# when they are installed; the core logistic/Cox/Poisson workflow depends on
# base R and (for Cox models) the survival package.

.r4vn_score_null <- function(x, y) if (is.null(x) || length(x) == 0L) y else x

.r4vn_score_bt <- function(x) {
  x <- as.character(x)
  paste0("`", gsub("`", "", x, fixed = TRUE), "`")
}

.r4vn_score_num <- function(x, digits = 3L) {
  ifelse(is.na(x), NA_character_, formatC(x, format = "f", digits = digits))
}

.r4vn_score_pct <- function(x, digits = 1L) {
  ifelse(is.na(x), NA_character_, paste0(formatC(100 * x, format = "f", digits = digits), "%"))
}

.r4vn_score_ci_text <- function(est, lo, hi, digits = 3L) {
  ifelse(
    is.na(est), NA_character_,
    paste0(.r4vn_score_num(est, digits), " (", .r4vn_score_num(lo, digits), "\u2013", .r4vn_score_num(hi, digits), ")")
  )
}

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

.r4vn_score_active_data <- function() {
  # Try common internal helpers first. This deliberately avoids requiring a
  # specific active-data implementation so the file can coexist with older
  # and newer R4VN releases.
  for (fn in c(".r4vn_get_active", ".r4vn_get_active_data", "active_data", "get_active_data", ".r4vn_active_data", "r4vn_active_data")) {
    if (exists(fn, mode = "function", inherits = TRUE)) {
      z <- try(get(fn, mode = "function", inherits = TRUE)(), silent = TRUE)
      if (!inherits(z, "try-error") && is.data.frame(z)) return(z)
    }
  }
  for (op in c("R4VN.active_data", "r4vn.active_data")) {
    z <- getOption(op, NULL)
    if (is.data.frame(z)) return(z)
    if (is.character(z) && length(z) == 1L && exists(z, envir = .GlobalEnv, inherits = FALSE)) {
      d <- get(z, envir = .GlobalEnv, inherits = FALSE)
      if (is.data.frame(d)) return(d)
    }
  }
  stop("No data supplied and no active R4VN data frame could be found. Supply data= or call usedf() first.", call. = FALSE)
}

.r4vn_score_extract_model <- function(x) {
  if (inherits(x, "glm") || inherits(x, "coxph")) return(x)
  if (is.list(x)) {
    candidates <- list(
      if (!is.null(x$raw) && is.list(x$raw)) x$raw$model else NULL,
      x$model,
      x$multi_model
    )
    for (z in candidates) if (inherits(z, "glm") || inherits(z, "coxph")) return(z)
  }
  NULL
}

.r4vn_score_resolve_name <- function(expr, env) {
  if (is.symbol(expr)) return(as.character(expr))
  val <- try(eval(expr, env), silent = TRUE)
  if (!inherits(val, "try-error") && is.character(val) && length(val) == 1L) return(val)
  if (is.character(expr) && length(expr) == 1L) return(expr)
  deparse(expr, width.cutoff = 500L)
}

.r4vn_score_parse_token <- function(txt, default_type = NA_character_) {
  txt <- as.character(txt)[1L]
  txt <- trimws(txt)
  ref <- NA_integer_
  type <- default_type
  if (grepl("^b[1-9][0-9]*\\.", txt)) {
    ref <- as.integer(sub("^b([1-9][0-9]*)\\..*$", "\\1", txt))
    txt <- sub("^b[1-9][0-9]*\\.", "", txt)
    type <- "categorical"
  } else if (grepl("^[cqf]\\.", txt)) {
    txt <- sub("^[cqf]\\.", "", txt)
    type <- "continuous"
  }
  data.frame(variable = txt, declared_type = type, reference_index = ref,
             stringsAsFactors = FALSE)
}

.r4vn_score_predictor_spec <- function(expr, env) {
  if (missing(expr) || is.null(expr) || identical(expr, quote(NULL))) {
    return(data.frame(variable = character(), declared_type = character(),
                      reference_index = integer(), stringsAsFactors = FALSE))
  }
  val <- try(eval(expr, env), silent = TRUE)
  if (!inherits(val, "try-error")) {
    # Native R4VN vars() object/data frame: preserve both declaration type and
    # requested categorical reference index rather than keeping only names.
    if (is.list(val) && !is.null(val$variable) && is.character(val$variable)) {
      type_source <- NULL
      for (nm in c("type", "kind", "summary", "mode")) {
        if (!is.null(val[[nm]]) && length(val[[nm]]) == length(val$variable)) {
          type_source <- as.character(val[[nm]])
          break
        }
      }
      typ0 <- if (is.null(type_source)) rep(NA_character_, length(val$variable)) else tolower(type_source)
      typ <- ifelse(grepl("categor|factor|binary|nominal|ordinal", typ0), "categorical",
                    ifelse(grepl("mean|median|full|continuous|numeric|quant", typ0), "continuous", NA_character_))
      ref_source <- NULL
      for (nm in c("reference_index", "ref_index", "reference", "ref")) {
        if (!is.null(val[[nm]]) && length(val[[nm]]) == length(val$variable)) {
          ref_source <- val[[nm]]
          break
        }
      }
      ref <- if (is.null(ref_source)) rep(NA_integer_, length(val$variable)) else suppressWarnings(as.integer(ref_source))
      z <- data.frame(variable = as.character(val$variable), declared_type = typ,
                      reference_index = ref, stringsAsFactors = FALSE)
      return(z[!duplicated(z$variable), , drop = FALSE])
    }
    if (is.character(val)) {
      z <- do.call(rbind, lapply(val, .r4vn_score_parse_token, default_type = NA_character_))
      return(z[!duplicated(z$variable), , drop = FALSE])
    }
    if (is.name(val)) return(.r4vn_score_parse_token(as.character(val), NA_character_))
  }
  if (is.call(expr)) {
    fn <- as.character(expr[[1L]])
    if (fn %in% c("vars", "c")) {
      default_type <- if (fn == "vars") "categorical" else NA_character_
      items <- as.list(expr)[-1L]
      z <- do.call(rbind, lapply(items, function(x) {
        .r4vn_score_parse_token(deparse(x, width.cutoff = 500L), default_type)
      }))
      return(z[!duplicated(z$variable), , drop = FALSE])
    }
  }
  txt <- deparse(expr, width.cutoff = 500L)
  txt <- gsub("^c\\(|^vars\\(|\\)$", "", txt)
  vals <- trimws(strsplit(txt, ",", fixed = TRUE)[[1L]])
  z <- do.call(rbind, lapply(vals, .r4vn_score_parse_token, default_type = NA_character_))
  z[!duplicated(z$variable), , drop = FALSE]
}

.r4vn_score_predictor_names <- function(expr, env) {
  .r4vn_score_predictor_spec(expr, env)$variable
}

.r4vn_score_prepare_declared <- function(data, spec) {
  out <- as.data.frame(data)
  if (is.null(spec) || !nrow(spec)) return(out)
  for (i in seq_len(nrow(spec))) {
    v <- spec$variable[[i]]
    if (!v %in% names(out)) next
    typ <- spec$declared_type[[i]]
    ref <- spec$reference_index[[i]]
    x <- out[[v]]
    if (identical(typ, "categorical")) {
      if (is.factor(x)) {
        f <- droplevels(x)
      } else if (is.character(x)) {
        lev <- unique(x[!is.na(x)])
        f <- factor(x, levels = lev)
      } else if (is.logical(x)) {
        f <- factor(ifelse(is.na(x), NA, ifelse(x, "Yes", "No")), levels = c("No", "Yes"))
      } else {
        lev <- sort(unique(x[!is.na(x)]))
        value_labels <- attr(x, "labels", exact = TRUE)
        if (!is.null(value_labels) && length(value_labels) && !is.null(names(value_labels)) &&
            all(nzchar(names(value_labels))) && all(lev %in% as.numeric(value_labels))) {
          labs <- names(value_labels)[match(lev, as.numeric(value_labels))]
        } else {
          labs <- format(lev, trim = TRUE, scientific = FALSE)
        }
        f <- factor(x, levels = lev, labels = labs)
      }
      if (is.finite(ref)) {
        if (ref < 1L || ref > nlevels(f))
          stop("Reference index b", ref, ". for predictor '", v,
               "' exceeds its ", nlevels(f), " observed level(s).", call. = FALSE)
        lv <- levels(f)
        f <- factor(as.character(f), levels = c(lv[[ref]], lv[-ref]))
      }
      out[[v]] <- f
    } else if (identical(typ, "continuous")) {
      if (!is.numeric(x))
        stop("Predictor '", v, "' was declared continuous (c./q./f.) but is not numeric.", call. = FALSE)
      out[[v]] <- as.numeric(x)
    }
  }
  out
}

.r4vn_score_capture_model_input <- function(data, predictors) {
  out <- list()
  for (v in predictors) {
    x <- data[[v]]
    if (is.factor(x)) out[[v]] <- list(type = "factor", levels = levels(x))
    else if (is.numeric(x)) out[[v]] <- list(type = "numeric")
    else if (is.logical(x)) out[[v]] <- list(type = "logical")
    else out[[v]] <- list(type = "character")
  }
  out
}

.r4vn_score_apply_model_input <- function(data, spec) {
  out <- as.data.frame(data)
  if (is.null(spec) || !length(spec)) return(out)
  for (v in names(spec)) {
    if (!v %in% names(out)) stop("newdata is missing original-model predictor '", v, "'.", call. = FALSE)
    z <- spec[[v]]
    if (identical(z$type, "factor")) {
      out[[v]] <- factor(as.character(out[[v]]), levels = z$levels)
    } else if (identical(z$type, "numeric")) {
      out[[v]] <- as.numeric(out[[v]])
    } else if (identical(z$type, "logical")) {
      out[[v]] <- as.logical(out[[v]])
    } else {
      out[[v]] <- as.character(out[[v]])
    }
  }
  out
}

.r4vn_score_rebuild_predictor_tokens <- function(spec) {
  if (is.null(spec) || !nrow(spec)) return(character())
  vapply(seq_len(nrow(spec)), function(i) {
    v <- spec$variable[[i]]
    typ <- spec$declared_type[[i]]
    ref <- spec$reference_index[[i]]
    if (identical(typ, "categorical")) {
      # Rebuild data are already prepared with the requested reference first.
      # b1. preserves categorical treatment without re-applying the original
      # positional reference index to an already releveled factor.
      paste0("b1.", v)
    } else if (identical(typ, "continuous")) paste0("c.", v) else v
  }, character(1L))
}

.r4vn_score_binary <- function(y, event = NULL, name = "outcome") {
  old <- y
  if (is.logical(y)) y <- as.integer(y)
  if (is.factor(y) || is.character(y)) {
    yy <- as.character(y)
    lev <- if (is.factor(y)) levels(droplevels(y)) else unique(yy[!is.na(yy)])
    if (length(lev) != 2L) stop(name, " must have exactly two non-missing levels.", call. = FALSE)
    if (is.null(event)) event <- lev[[2L]]
    if (!as.character(event) %in% lev) stop("event='", event, "' was not found in ", name, ".", call. = FALSE)
    y <- ifelse(is.na(yy), NA_integer_, as.integer(yy == as.character(event)))
    return(list(y = y, event = as.character(event), levels = lev, original = old))
  }
  u <- sort(unique(y[!is.na(y)]))
  if (length(u) != 2L) stop(name, " must be binary for this analysis.", call. = FALSE)
  if (is.null(event)) {
    if (all(u %in% c(0, 1))) event <- 1 else event <- u[[2L]]
  }
  if (!event %in% u) stop("event= was not found in ", name, ".", call. = FALSE)
  list(y = ifelse(is.na(y), NA_integer_, as.integer(y == event)), event = event, levels = u, original = old)
}

.r4vn_score_is_binary <- function(y) {
  length(unique(y[!is.na(y)])) == 2L
}

.r4vn_score_formula <- function(family, predictors, outcome = ".r4vn_score_y", time = ".r4vn_score_time") {
  rhs <- if (length(predictors)) paste(.r4vn_score_bt(predictors), collapse = " + ") else "1"
  if (family == "cox") {
    stats::as.formula(paste0("survival::Surv(", time, ", ", outcome, ") ~ ", rhs))
  } else {
    stats::as.formula(paste0(outcome, " ~ ", rhs))
  }
}

.r4vn_score_fit <- function(data, family, predictors) {
  f <- .r4vn_score_formula(family, predictors)
  if (family == "logistic") {
    stats::glm(f, data = data, family = stats::binomial(link = "logit"), x = TRUE, y = TRUE, model = TRUE)
  } else if (family == "poisson") {
    stats::glm(f, data = data, family = stats::poisson(link = "log"), x = TRUE, y = TRUE, model = TRUE)
  } else if (family == "cox") {
    if (!requireNamespace("survival", quietly = TRUE)) stop("Package 'survival' is required for Cox scorecards.", call. = FALSE)
    survival::coxph(f, data = data, x = TRUE, y = TRUE, model = TRUE, ties = "efron")
  } else {
    stop("Unsupported family: ", family, call. = FALSE)
  }
}

.r4vn_score_lrt <- function(smaller, larger, family) {
  out <- try({
    a <- if (family == "cox") stats::anova(smaller, larger, test = "Chisq") else stats::anova(smaller, larger, test = "Chisq")
    pcol <- grep("Pr\\(", colnames(a), value = TRUE)
    if (!length(pcol)) NA_real_ else as.numeric(a[nrow(a), pcol[[1L]]])
  }, silent = TRUE)
  if (inherits(out, "try-error") || !is.finite(out)) NA_real_ else out
}

.r4vn_score_term_p <- function(data, family, current, term) {
  if (!term %in% current) return(NA_real_)
  full <- .r4vn_score_fit(data, family, current)
  reduced <- .r4vn_score_fit(data, family, setdiff(current, term))
  .r4vn_score_lrt(reduced, full, family)
}

.r4vn_score_coef_change <- function(full, reduced) {
  b1 <- stats::coef(full)
  b0 <- stats::coef(reduced)
  nm <- intersect(names(b1), names(b0))
  nm <- setdiff(nm, "(Intercept)")
  if (!length(nm)) return(0)
  den <- pmax(abs(b1[nm]), 1e-6)
  max(abs(b0[nm] - b1[nm]) / den, na.rm = TRUE)
}

.r4vn_score_select_purposeful <- function(data, family, predictors, force = NULL,
                                           entry = 0.25, stay = 0.10, confound = 0.15) {
  force <- intersect(.r4vn_score_null(force, character()), predictors)
  null <- .r4vn_score_fit(data, family, character())
  up <- vapply(predictors, function(v) {
    fit <- try(.r4vn_score_fit(data, family, v), silent = TRUE)
    if (inherits(fit, "try-error")) return(NA_real_)
    .r4vn_score_lrt(null, fit, family)
  }, numeric(1L))
  current <- unique(c(force, names(up)[is.na(up) | up < entry]))
  if (!length(current)) current <- force
  if (!length(current)) current <- predictors[[which.min(replace(up, is.na(up), Inf))]]

  protected <- force
  repeat {
    removable <- setdiff(current, protected)
    if (!length(removable)) break
    ps <- vapply(removable, function(v) .r4vn_score_term_p(data, family, current, v), numeric(1L))
    if (all(is.na(ps))) break
    worst <- removable[[which.max(replace(ps, is.na(ps), -Inf))]]
    pworst <- ps[[worst]]
    if (!is.finite(pworst) || pworst <= stay) break
    full <- .r4vn_score_fit(data, family, current)
    reduced_terms <- setdiff(current, worst)
    reduced <- .r4vn_score_fit(data, family, reduced_terms)
    chg <- .r4vn_score_coef_change(full, reduced)
    if (is.finite(chg) && chg > confound) {
      protected <- unique(c(protected, worst))
    } else {
      current <- reduced_terms
    }
  }

  # Re-assess variables excluded at screening, one at a time.
  excluded <- setdiff(predictors, current)
  if (length(excluded)) {
    base <- .r4vn_score_fit(data, family, current)
    for (v in excluded) {
      add <- try(.r4vn_score_fit(data, family, c(current, v)), silent = TRUE)
      if (!inherits(add, "try-error")) {
        p <- .r4vn_score_lrt(base, add, family)
        if (is.finite(p) && p < stay) {
          current <- c(current, v)
          base <- add
        }
      }
    }
  }
  unique(c(force, current))
}

.r4vn_score_select_lasso <- function(data, family, predictors, force = NULL, seed = NULL) {
  if (!requireNamespace("glmnet", quietly = TRUE)) {
    stop("select='lasso' requires package 'glmnet'. Install it or choose another selection method.", call. = FALSE)
  }
  f <- .r4vn_score_formula(family, predictors)
  mm <- stats::model.matrix(f, data = data)
  if ("(Intercept)" %in% colnames(mm)) mm <- mm[, colnames(mm) != "(Intercept)", drop = FALSE]
  if (!ncol(mm)) return(unique(force))
  yy <- if (family == "cox") survival::Surv(data$.r4vn_score_time, data$.r4vn_score_y) else data$.r4vn_score_y
  fam <- if (family == "logistic") "binomial" else if (family == "cox") "cox" else "poisson"
  .r4vn_set_seed_if(seed)
  cv <- glmnet::cv.glmnet(mm, yy, family = fam, alpha = 1)
  cc <- as.matrix(stats::coef(cv, s = "lambda.1se"))
  nz <- rownames(cc)[as.numeric(cc[, 1L]) != 0]
  nz <- setdiff(nz, "(Intercept)")
  if (!length(nz)) {
    cc <- as.matrix(stats::coef(cv, s = "lambda.min"))
    nz <- rownames(cc)[as.numeric(cc[, 1L]) != 0]
    nz <- setdiff(nz, "(Intercept)")
  }
  # Map selected dummy columns back to original predictor terms.
  assign <- attr(stats::model.matrix(f, data = data), "assign")
  full_mm_names <- colnames(stats::model.matrix(f, data = data))
  term_labels <- attr(stats::terms(f), "term.labels")
  selected <- character()
  for (z in nz) {
    j <- match(z, full_mm_names)
    if (is.finite(j) && !is.na(j) && assign[[j]] > 0L) selected <- c(selected, term_labels[[assign[[j]]]])
  }
  selected <- unique(gsub("`", "", selected, fixed = TRUE))
  selected <- intersect(selected, predictors)
  unique(c(intersect(force, predictors), selected))
}

.r4vn_score_select <- function(data, family, predictors, method = "none", force = NULL,
                               entry = 0.25, stay = 0.10, confound = 0.15, seed = NULL) {
  method <- match.arg(method, c("none", "full", "backward", "forward", "purposeful", "lasso"))
  force <- intersect(.r4vn_score_null(force, character()), predictors)
  if (method %in% c("none", "full")) return(unique(c(force, predictors)))
  if (method == "purposeful") {
    return(.r4vn_score_select_purposeful(data, family, predictors, force, entry, stay, confound))
  }
  if (method == "lasso") return(.r4vn_score_select_lasso(data, family, predictors, force, seed))

  full <- .r4vn_score_fit(data, family, predictors)
  if (method == "backward") {
    st <- stats::step(full, direction = "backward", trace = 0)
    trm <- attr(stats::terms(st), "term.labels")
    trm <- gsub("`", "", trm, fixed = TRUE)
    return(unique(c(force, intersect(trm, predictors))))
  }

  # Forward AIC selection is implemented explicitly rather than through
  # stats::step().  A fitted glm stores the original call as `data = data`;
  # when step() later evaluates that call from its own frame, the local
  # data-frame binding may no longer be available (notably in examples and
  # package checks).  Re-fitting candidates from the data object passed here
  # avoids that evaluation-environment dependency and also keeps forced terms
  # in every candidate model.
  current <- unique(force)
  current_fit <- .r4vn_score_fit(data, family, current)
  current_aic <- suppressWarnings(try(stats::AIC(current_fit), silent = TRUE))
  current_aic <- if (inherits(current_aic, "try-error") || !length(current_aic)) {
    Inf
  } else {
    as.numeric(current_aic[[1L]])
  }

  remaining <- setdiff(predictors, current)
  while (length(remaining)) {
    fits <- lapply(remaining, function(v) {
      try(.r4vn_score_fit(data, family, c(current, v)), silent = TRUE)
    })
    aic <- vapply(fits, function(fit) {
      if (inherits(fit, "try-error")) return(Inf)
      z <- suppressWarnings(try(stats::AIC(fit), silent = TRUE))
      if (inherits(z, "try-error") || !length(z)) return(Inf)
      z <- as.numeric(z[[1L]])
      if (is.finite(z)) z else Inf
    }, numeric(1L))

    if (!length(aic) || all(!is.finite(aic))) break
    j <- which.min(aic)
    # Match step()'s practical rule: add a term only when AIC is genuinely
    # improved, allowing a tiny numerical tolerance for equal models.
    if (!is.finite(aic[[j]]) || aic[[j]] >= current_aic - 1e-7) break

    current <- c(current, remaining[[j]])
    current_fit <- fits[[j]]
    current_aic <- aic[[j]]
    remaining <- setdiff(predictors, current)
  }

  unique(current)
}

.r4vn_score_nice_step <- function(x) {
  x <- abs(x)
  x <- x[is.finite(x) & x > 0]
  if (!length(x)) return(1)
  target <- stats::median(x)
  power <- 10^floor(log10(target))
  z <- target / power
  mult <- if (z < 1.5) 1 else if (z < 3.5) 2 else if (z < 7.5) 5 else 10
  mult * power
}

.r4vn_score_easy_cuts <- function(x, bins = 4L) {
  x <- x[is.finite(x)]
  if (length(unique(x)) < 4L) return(numeric())
  probs <- seq(0, 1, length.out = bins + 1L)[-c(1L, bins + 1L)]
  q <- as.numeric(stats::quantile(x, probs = probs, na.rm = TRUE, names = FALSE, type = 2))
  rng <- diff(range(x, na.rm = TRUE))
  step <- .r4vn_score_nice_step(rng / (bins * 2))
  if (!is.finite(step) || step <= 0) step <- 1
  nice <- unique(round(q / step) * step)
  nice <- nice[nice > min(x, na.rm = TRUE) & nice < max(x, na.rm = TRUE)]
  if (length(nice) < 1L) nice <- unique(q[q > min(x) & q < max(x)])
  nice
}

.r4vn_score_cut_labels <- function(x, cuts) {
  cuts <- sort(unique(as.numeric(cuts)))
  integer_like <- all(abs(x[is.finite(x)] - round(x[is.finite(x)])) < 1e-8) &&
    all(abs(cuts - round(cuts)) < 1e-8)
  f <- function(z) format(z, trim = TRUE, scientific = FALSE)
  if (!length(cuts)) return(character())
  labs <- character(length(cuts) + 1L)
  labs[[1L]] <- paste0("<", f(cuts[[1L]]))
  if (length(cuts) > 1L) {
    for (i in 2:length(cuts)) {
      lo <- cuts[[i - 1L]]
      hi <- cuts[[i]]
      if (integer_like) labs[[i]] <- paste0(f(lo), "\u2013", f(hi - 1)) else labs[[i]] <- paste0(f(lo), "\u2013<", f(hi))
    }
  }
  labs[[length(labs)]] <- paste0("\u2265", f(cuts[[length(cuts)]]))
  labs
}

.r4vn_score_transform_fit <- function(data, predictors, cuts = NULL,
                                      continuous = "easy", bins = 4L) {
  continuous <- match.arg(continuous, c("easy", "auto", "quantile", "keep"))
  out <- data
  spec <- list()
  for (v in predictors) {
    x <- data[[v]]
    manual <- is.list(cuts) && !is.null(cuts[[v]])
    if (manual) {
      br <- sort(unique(as.numeric(cuts[[v]])))
      labs <- .r4vn_score_cut_labels(x, br)
      out[[v]] <- cut(x, breaks = c(-Inf, br, Inf), labels = labs, right = FALSE, include.lowest = TRUE)
      spec[[v]] <- list(type = "cut", breaks = br, labels = labs)
      next
    }
    if (is.character(x)) x <- factor(x)
    if (is.logical(x)) x <- factor(ifelse(is.na(x), NA, ifelse(x, "Yes", "No")), levels = c("No", "Yes"))
    value_labels <- attr(x, "labels", exact = TRUE)
    if (is.numeric(x) && !is.null(value_labels) && length(value_labels) &&
        !is.null(names(value_labels)) && all(nzchar(names(value_labels)))) {
      obs <- sort(unique(as.numeric(x[!is.na(x)])))
      vv <- as.numeric(value_labels)
      if (length(obs) && all(obs %in% vv)) {
        ll <- names(value_labels)[match(obs, vv)]
        out[[v]] <- factor(as.numeric(x), levels = obs, labels = ll)
        spec[[v]] <- list(type = "factor_numeric", levels = obs, labels = ll)
        next
      }
    }
    if (is.factor(x)) {
      out[[v]] <- droplevels(x)
      spec[[v]] <- list(type = "factor", levels = levels(out[[v]]))
      next
    }
    if (is.numeric(x)) {
      u <- sort(unique(x[!is.na(x)]))
      if (length(u) <= 10L && continuous != "keep") {
        out[[v]] <- factor(x, levels = u, labels = format(u, trim = TRUE, scientific = FALSE))
        spec[[v]] <- list(type = "factor_numeric", levels = u, labels = levels(out[[v]]))
      } else if (continuous == "keep") {
        out[[v]] <- x
        spec[[v]] <- list(type = "numeric", range = range(x, na.rm = TRUE))
      } else {
        br <- if (continuous == "quantile") {
          unique(as.numeric(stats::quantile(x, probs = seq(0, 1, length.out = bins + 1L)[-c(1L, bins + 1L)],
                                             na.rm = TRUE, names = FALSE, type = 2)))
        } else .r4vn_score_easy_cuts(x, bins)
        br <- br[br > min(x, na.rm = TRUE) & br < max(x, na.rm = TRUE)]
        if (!length(br)) {
          out[[v]] <- x
          spec[[v]] <- list(type = "numeric", range = range(x, na.rm = TRUE))
        } else {
          labs <- .r4vn_score_cut_labels(x, br)
          out[[v]] <- cut(x, breaks = c(-Inf, br, Inf), labels = labs, right = FALSE, include.lowest = TRUE)
          spec[[v]] <- list(type = "cut", breaks = br, labels = labs)
        }
      }
      next
    }
    stop("Predictor '", v, "' has an unsupported class for score construction.", call. = FALSE)
  }
  list(data = out, spec = spec)
}

.r4vn_score_transform_apply <- function(data, spec) {
  out <- data
  for (v in names(spec)) {
    if (!v %in% names(out)) stop("newdata is missing predictor '", v, "'.", call. = FALSE)
    s <- spec[[v]]
    x <- out[[v]]
    if (s$type == "cut") {
      out[[v]] <- cut(as.numeric(x), breaks = c(-Inf, s$breaks, Inf), labels = s$labels,
                      right = FALSE, include.lowest = TRUE)
    } else if (s$type == "factor") {
      out[[v]] <- factor(as.character(x), levels = s$levels)
    } else if (s$type == "factor_numeric") {
      out[[v]] <- factor(as.numeric(x), levels = s$levels, labels = s$labels)
    } else if (s$type == "numeric") {
      out[[v]] <- as.numeric(x)
    }
  }
  out
}

.r4vn_score_reference_row <- function(data, predictors) {
  z <- lapply(predictors, function(v) {
    x <- data[[v]]
    if (is.factor(x)) factor(levels(x)[1L], levels = levels(x)) else if (is.numeric(x)) stats::median(x, na.rm = TRUE) else x[[which(!is.na(x))[1L]]]
  })
  names(z) <- predictors
  as.data.frame(z, check.names = FALSE)
}

.r4vn_score_effect_measure <- function(family, y = NULL) {
  if (identical(family, "logistic")) return("OR")
  if (identical(family, "cox")) return("HR")
  if (identical(family, "poisson")) {
    if (!is.null(y) && .r4vn_score_is_binary(y)) return("RR")
    return("IRR")
  }
  "Effect ratio"
}

.r4vn_score_scoreref_value <- function(scoreref, predictor) {
  if (is.null(scoreref)) return(NULL)
  if (is.list(scoreref)) {
    if (is.null(names(scoreref)) || !predictor %in% names(scoreref)) return(NULL)
    z <- scoreref[[predictor]]
    if (is.null(z) || !length(z) || is.na(z[[1L]])) return(NULL)
    return(as.character(z[[1L]]))
  }
  if (is.character(scoreref) && !is.null(names(scoreref))) {
    if (!predictor %in% names(scoreref)) return(NULL)
    z <- scoreref[[predictor]]
    if (is.na(z) || !nzchar(z)) return(NULL)
    return(as.character(z))
  }
  NULL
}

.r4vn_score_validate_scoreref <- function(scoreref) {
  if (is.null(scoreref)) return(invisible(TRUE))
  if (is.character(scoreref) && length(scoreref) == 1L && is.null(names(scoreref))) {
    if (!tolower(scoreref) %in% c("lowest", "model"))
      stop("`scoreref` must be 'lowest', 'model', or a named list/vector of score-reference categories.", call. = FALSE)
    return(invisible(TRUE))
  }
  if (is.list(scoreref) || (is.character(scoreref) && !is.null(names(scoreref)))) {
    if (is.null(names(scoreref)) || any(!nzchar(names(scoreref))))
      stop("A custom `scoreref` must be named by predictor, e.g. list(sex='Female').", call. = FALSE)
    return(invisible(TRUE))
  }
  stop("`scoreref` must be 'lowest', 'model', or a named list/vector of score-reference categories.", call. = FALSE)
}

.r4vn_score_predict_lp <- function(fit, newdata = NULL) {
  # glm models use type = "link". survival::coxph uses type = "lp".
  # reference = "zero" keeps Cox predictions on the X beta scale so that
  # category contributions use one common origin across prediction calls.
  if (inherits(fit, "coxph")) {
    args <- list(object = fit, type = "lp", reference = "zero")
    if (!is.null(newdata)) args$newdata <- newdata
    return(as.numeric(do.call(stats::predict, args)))
  }
  args <- list(object = fit, type = "link")
  if (!is.null(newdata)) args$newdata <- newdata
  as.numeric(do.call(stats::predict, args))
}

.r4vn_score_effect_dictionary <- function(fit, data, predictors, labels,
                                          family, riskonly = TRUE,
                                          scoreref = "lowest") {
  .r4vn_score_validate_scoreref(scoreref)
  ref <- .r4vn_score_reference_row(data, predictors)
  lp_ref <- .r4vn_score_predict_lp(fit, newdata = ref)
  rows <- list()
  numeric_present <- character()
  tol <- 1e-10

  global_mode <- if (is.character(scoreref) && length(scoreref) == 1L && is.null(names(scoreref))) {
    tolower(scoreref)
  } else {
    "custom"
  }

  for (v in predictors) {
    x <- data[[v]]
    if (is.factor(x)) {
      levs <- levels(x)
      raw_eff <- vapply(levs, function(lev) {
        nd <- ref
        nd[[v]] <- factor(lev, levels = levs)
        .r4vn_score_predict_lp(fit, newdata = nd) - lp_ref
      }, numeric(1L))
      names(raw_eff) <- levs

      model_ref_index <- 1L
      requested <- .r4vn_score_scoreref_value(scoreref, v)
      if (!is.null(requested)) {
        score_ref_index <- match(requested, levs)
        if (is.na(score_ref_index)) {
          stop("Scoring reference '", requested, "' for predictor '", v,
               "' was not found. Available categories: ", paste(levs, collapse = ", "), ".",
               call. = FALSE)
        }
      } else if (identical(global_mode, "model")) {
        score_ref_index <- model_ref_index
      } else {
        mn <- min(raw_eff, na.rm = TRUE)
        low <- which(abs(raw_eff - mn) <= tol)
        score_ref_index <- if (model_ref_index %in% low) model_ref_index else low[[1L]]
      }

      if (isTRUE(riskonly)) {
        mn <- min(raw_eff, na.rm = TRUE)
        if (raw_eff[[score_ref_index]] > mn + tol) {
          stop(
            "`riskonly=TRUE` requires the scoring reference for predictor '", v,
            "' to be a lowest-risk category. Requested '", levs[[score_ref_index]],
            "' is not the lowest-risk category. Use `scoreref='lowest'` (recommended), ",
            "choose a lowest-risk category explicitly, or set `riskonly=FALSE` if negative points are intended.",
            call. = FALSE
          )
        }
        # Rebase each predictor to its minimum modeled contribution. This is a
        # reparameterization: the predictor contrast ratios become >= 1 while
        # the full model predictions/ranking remain unchanged.
        score_eff <- raw_eff - mn
      } else {
        score_eff <- raw_eff - raw_eff[[score_ref_index]]
      }

      # Numerical noise around zero should never create tiny negative points.
      score_eff[abs(score_eff) <= tol] <- 0
      if (isTRUE(riskonly) && any(score_eff < -tol, na.rm = TRUE)) {
        stop("Internal score orientation failed for predictor '", v, "'.", call. = FALSE)
      }

      rows[[v]] <- data.frame(
        predictor = v,
        predictor_label = labels[[v]],
        category = levs,
        model_effect = as.numeric(raw_eff),
        model_ratio = exp(as.numeric(raw_eff)),
        effect = as.numeric(score_eff),
        scoring_ratio = exp(as.numeric(score_eff)),
        model_reference = seq_along(levs) == model_ref_index,
        scoring_reference = seq_along(levs) == score_ref_index,
        zero_point_category = abs(score_eff) <= tol,
        protective_vs_model_reference = raw_eff < -tol,
        scoring_rule = "Category",
        stringsAsFactors = FALSE
      )
    } else if (is.numeric(x)) {
      numeric_present <- c(numeric_present, v)
      b <- stats::coef(fit)
      bn <- names(b)
      j <- match(v, gsub("`", "", bn, fixed = TRUE))
      beta <- if (!is.na(j)) as.numeric(b[[j]]) else NA_real_
      score_beta <- if (isTRUE(riskonly)) abs(beta) else beta
      rule <- if (!is.finite(beta)) {
        "Per 1 unit"
      } else if (isTRUE(riskonly) && beta < 0) {
        "Per 1-unit decrease"
      } else {
        "Per 1-unit increase"
      }
      rows[[v]] <- data.frame(
        predictor = v,
        predictor_label = labels[[v]],
        category = "Continuous",
        model_effect = beta,
        model_ratio = exp(beta),
        effect = score_beta,
        scoring_ratio = exp(score_beta),
        model_reference = FALSE,
        scoring_reference = FALSE,
        zero_point_category = FALSE,
        protective_vs_model_reference = is.finite(beta) && beta < 0,
        scoring_rule = rule,
        stringsAsFactors = FALSE
      )
    }
  }
  tab <- do.call(rbind, rows)
  if (!is.null(tab) && nrow(tab)) rownames(tab) <- NULL
  list(
    table = tab,
    numeric = numeric_present,
    effect_measure = .r4vn_score_effect_measure(family, data$.r4vn_score_y),
    riskonly = isTRUE(riskonly),
    scoreref = scoreref
  )
}

.r4vn_score_auc <- function(y, marker) {
  ok <- is.finite(marker) & !is.na(y)
  y <- y[ok]; marker <- marker[ok]
  if (length(unique(y)) != 2L) return(NA_real_)
  n1 <- sum(y == 1); n0 <- sum(y == 0)
  if (!n1 || !n0) return(NA_real_)
  r <- rank(marker, ties.method = "average")
  (sum(r[y == 1]) - n1 * (n1 + 1) / 2) / (n1 * n0)
}

.r4vn_score_auc_ci <- function(y, marker, conf.level = 0.95, B = 300L, seed = NULL) {
  est <- .r4vn_score_auc(y, marker)
  if (requireNamespace("pROC", quietly = TRUE) && is.finite(est)) {
    rr <- try(pROC::roc(response = y, predictor = marker, quiet = TRUE, direction = "<"), silent = TRUE)
    if (!inherits(rr, "try-error")) {
      ci <- try(as.numeric(pROC::ci.auc(rr, conf.level = conf.level)), silent = TRUE)
      if (!inherits(ci, "try-error") && length(ci) == 3L) return(c(est = est, low = ci[[1L]], high = ci[[3L]]))
    }
  }
  .r4vn_set_seed_if(seed)
  n <- length(y)
  bb <- replicate(B, {
    ii <- sample.int(n, n, replace = TRUE)
    .r4vn_score_auc(y[ii], marker[ii])
  })
  a <- (1 - conf.level) / 2
  bb <- bb[is.finite(bb)]
  if (!length(bb)) return(c(est = est, low = NA_real_, high = NA_real_))
  ci <- stats::quantile(bb, c(a, 1 - a), na.rm = TRUE, names = FALSE)
  c(est = est, low = ci[[1L]], high = ci[[2L]])
}

.r4vn_score_cindex <- function(time, event, marker, conf.level = 0.95) {
  if (!requireNamespace("survival", quietly = TRUE)) return(c(est = NA, low = NA, high = NA))
  cc <- try(survival::concordance(survival::Surv(time, event) ~ marker, reverse = TRUE), silent = TRUE)
  if (inherits(cc, "try-error")) return(c(est = NA, low = NA, high = NA))
  est <- as.numeric(cc$concordance)
  se <- sqrt(as.numeric(cc$var))
  z <- stats::qnorm(1 - (1 - conf.level) / 2)
  c(est = est, low = max(0, est - z * se), high = min(1, est + z * se))
}

.r4vn_score_perf_metric <- function(family, data, marker) {
  if (family == "logistic" || (family == "poisson" && .r4vn_score_is_binary(data$.r4vn_score_y))) {
    .r4vn_score_auc(data$.r4vn_score_y, marker)
  } else if (family == "cox") {
    .r4vn_score_cindex(data$.r4vn_score_time, data$.r4vn_score_y, marker)[["est"]]
  } else {
    suppressWarnings(stats::cor(marker, data$.r4vn_score_y, method = "spearman", use = "complete.obs"))
  }
}

.r4vn_score_apply_dictionary <- function(data, dictionary, point_col = "clinical_point") {
  s <- rep(0, nrow(data))
  for (v in unique(dictionary$predictor)) {
    d <- dictionary[dictionary$predictor == v, , drop = FALSE]
    x <- as.character(data[[v]])
    mp <- setNames(d[[point_col]], d$category)
    z <- unname(mp[x])
    s <- s + z
  }
  s
}

.r4vn_score_normalize_category_key <- function(x) {
  x <- as.character(x)
  x <- trimws(x)
  x <- gsub("\u00a0", " ", x, fixed = TRUE)
  x <- gsub("\u2265", ">=", x, fixed = TRUE)
  x <- gsub("\u2264", "<=", x, fixed = TRUE)
  x <- gsub("[\u2010\u2011\u2012\u2013\u2014\u2212]", "-", x, perl = TRUE)
  x <- gsub("\\s*([<>]=?)\\s*", "\\1", x, perl = TRUE)
  x <- gsub("\\s*-\\s*", "-", x, perl = TRUE)
  x <- gsub("\\s+", " ", x, perl = TRUE)
  trimws(x)
}

.r4vn_score_match_categories <- function(expected, supplied, predictor) {
  expected <- as.character(expected)
  supplied <- as.character(supplied)
  expected_key <- .r4vn_score_normalize_category_key(expected)
  supplied_key <- .r4vn_score_normalize_category_key(supplied)

  dup <- unique(supplied_key[duplicated(supplied_key)])
  if (length(dup)) {
    stop(
      "Manual points for '", predictor,
      "' contain duplicated/ambiguous category names after normalization: ",
      paste(dup, collapse = ", "), ".",
      call. = FALSE
    )
  }

  idx <- match(expected_key, supplied_key)
  if (anyNA(idx)) {
    misscat <- expected[is.na(idx)]
    stop(
      "Manual points for '", predictor,
      "' are missing named category/categories: ",
      paste(misscat, collapse = ", "), ".",
      call. = FALSE
    )
  }
  idx
}

.r4vn_score_manual_points <- function(dict, points, riskonly = TRUE) {
  out <- dict
  out$clinical_point <- NA_real_
  if (is.data.frame(points)) {
    nms <- tolower(names(points))
    ip <- match("predictor", nms); ic <- match("category", nms); iv <- match("point", nms)
    if (anyNA(c(ip, ic, iv))) stop("Manual points data frame must contain Predictor, Category and Point columns.", call. = FALSE)

    point_predictor <- as.character(points[[ip]])
    for (v in unique(out$predictor)) {
      ii <- which(out$predictor == v)
      jj <- which(point_predictor == v)
      if (!length(jj)) stop("Manual points are missing predictor '", v, "'.", call. = FALSE)
      idx <- .r4vn_score_match_categories(
        expected = out$category[ii],
        supplied = points[[ic]][jj],
        predictor = v
      )
      out$clinical_point[ii] <- as.numeric(points[[iv]][jj][idx])
    }
  } else if (is.list(points)) {
    for (v in unique(out$predictor)) {
      if (is.null(points[[v]])) stop("Manual points are missing predictor '", v, "'.", call. = FALSE)
      d <- out[out$predictor == v, , drop = FALSE]
      p <- points[[v]]
      if (!is.null(names(p))) {
        idx <- .r4vn_score_match_categories(
          expected = d$category,
          supplied = names(p),
          predictor = v
        )
        p <- p[idx]
      }
      if (length(p) != nrow(d)) stop("Manual points for '", v, "' must have one value per category.", call. = FALSE)
      out$clinical_point[out$predictor == v] <- as.numeric(p)
    }
  } else stop("Manual points must be a named list or a data frame.", call. = FALSE)
  if (anyNA(out$clinical_point)) stop("Some manual points could not be matched to scorecard categories.", call. = FALSE)
  if (any(!is.finite(out$clinical_point)) || any(abs(out$clinical_point - round(out$clinical_point)) > 1e-8))
    stop("Manual clinical points must be finite integers.", call. = FALSE)
  out$clinical_point <- as.integer(round(out$clinical_point))

  if (isTRUE(riskonly)) {
    # Manual scores are also converted to add-only form predictor by predictor.
    # Subtracting the within-predictor minimum changes only the total-score
    # origin; it does not change subject ranking or any score contrast.
    for (v in unique(out$predictor)) {
      ii <- which(out$predictor == v)
      out$clinical_point[ii] <- out$clinical_point[ii] - min(out$clinical_point[ii], na.rm = TRUE)
    }
    if (any(out$clinical_point < 0L))
      stop("`riskonly=TRUE` could not convert manual points to a non-negative score.", call. = FALSE)
  }
  out
}

.r4vn_score_build_points <- function(dict, score_data, family, fit, pdo = 20,
                                     points = "auto", maxscore = "auto", tolerance = 0.01,
                                     riskonly = TRUE) {
  if (any(dict$category == "Continuous")) {
    stop("A complete bedside scorecard needs explicit score-ready categories for continuous predictors. Use continuous='easy' (default) or supply cuts=. continuous='keep' preserves the variable for modeling but cannot produce a complete Predictor\u2013Category\u2013Point table in this version.", call. = FALSE)
  }
  factor_pdo <- pdo / log(2)
  dict$model_point <- dict$effect * factor_pdo

  if (is.list(points) || is.data.frame(points)) {
    out <- .r4vn_score_manual_points(dict, points, riskonly = riskonly)
    return(list(dictionary = out, scale_B = NA_real_, method = "manual", warning = NULL))
  }
  points <- match.arg(points, c("auto", "clinical", "integer", "model", "pdo"))
  if (points %in% c("model", "pdo")) {
    out <- dict
    out$clinical_point <- round(out$model_point)
    if (isTRUE(riskonly)) {
      for (v in unique(out$predictor)) {
        ii <- which(out$predictor == v)
        out$clinical_point[ii] <- out$clinical_point[ii] - min(out$clinical_point[ii], na.rm = TRUE)
      }
    }
    return(list(dictionary = out, scale_B = log(2) / pdo, method = points, warning = NULL))
  }

  summax <- sum(vapply(split(dict$effect, dict$predictor), max, numeric(1L), na.rm = TRUE))
  if (!is.finite(summax) || summax <= 0) {
    out <- dict; out$clinical_point <- 0
    return(list(dictionary = out, scale_B = 1, method = "auto", warning = "All adjusted predictor contributions were zero."))
  }

  lp <- .r4vn_score_predict_lp(fit)
  target_metric <- .r4vn_score_perf_metric(family, score_data, lp)
  targets <- if (is.numeric(maxscore) && length(maxscore) == 1L) max(2L, round(maxscore)) else 5:30
  cand <- list(); k <- 0L
  for (tg in targets) {
    b0 <- summax / tg
    for (m in seq(0.75, 1.30, by = 0.025)) {
      B <- b0 * m
      pt <- round(dict$effect / B)
      if (all(pt == 0)) next
      dd <- dict; dd$clinical_point <- pt
      sc <- .r4vn_score_apply_dictionary(score_data, dd)
      mx <- sum(vapply(split(dd$clinical_point, dd$predictor), max, numeric(1L), na.rm = TRUE))
      mn <- sum(vapply(split(dd$clinical_point, dd$predictor), min, numeric(1L), na.rm = TRUE))
      if (is.numeric(maxscore) && mx - mn > maxscore) next
      met <- .r4vn_score_perf_metric(family, score_data, sc)
      if (!is.finite(met)) next
      k <- k + 1L
      cand[[k]] <- list(B = B, points = pt, score = sc, metric = met, range = mx - mn,
                        distinct = length(unique(pt)))
    }
  }
  if (!length(cand)) stop("Could not derive a non-zero clinical score. Consider a larger maxscore or manual points.", call. = FALSE)
  ctab <- data.frame(
    i = seq_along(cand),
    metric = vapply(cand, `[[`, numeric(1L), "metric"),
    range = vapply(cand, `[[`, numeric(1L), "range"),
    distinct = vapply(cand, `[[`, numeric(1L), "distinct")
  )
  ctab$loss <- target_metric - ctab$metric
  ok <- which(is.finite(ctab$loss) & ctab$loss <= tolerance)
  warn <- NULL
  if (length(ok)) {
    z <- ctab[ok, , drop = FALSE]
    z <- z[order(z$range, z$distinct, -z$metric), , drop = FALSE]
    best <- z$i[[1L]]
  } else {
    z <- ctab[order(-ctab$metric, ctab$range, ctab$distinct), , drop = FALSE]
    best <- z$i[[1L]]
    warn <- paste0("No candidate clinical score stayed within tolerance=", tolerance,
                   " of the score-ready model performance; the best-performing candidate was retained.")
  }
  out <- dict
  out$clinical_point <- cand[[best]]$points
  if (isTRUE(riskonly)) {
    for (v in unique(out$predictor)) {
      ii <- which(out$predictor == v)
      out$clinical_point[ii] <- out$clinical_point[ii] - min(out$clinical_point[ii], na.rm = TRUE)
    }
  }
  list(dictionary = out, scale_B = cand[[best]]$B, method = "auto",
       warning = warn, target_metric = target_metric, clinical_metric = cand[[best]]$metric,
       performance_loss = target_metric - cand[[best]]$metric)
}


.r4vn_score_attainable <- function(dict, point_col = "clinical_point") {
  vals <- split(dict[[point_col]], dict$predictor)
  sums <- 0
  for (v in vals) sums <- sort(unique(as.vector(outer(sums, unique(v), "+"))))
  sums[is.finite(sums)]
}

.r4vn_score_theoretical_range <- function(dict, point_col = "clinical_point") {
  sp <- split(dict[[point_col]], dict$predictor)
  c(min = sum(vapply(sp, min, numeric(1L), na.rm = TRUE)),
    max = sum(vapply(sp, max, numeric(1L), na.rm = TRUE)))
}

.r4vn_score_recal_binary <- function(y, score) {
  d <- data.frame(.y = y, .score = score)
  stats::glm(.y ~ .score, data = d, family = stats::binomial())
}

.r4vn_score_calibration <- function(y, p, conf.level = 0.95) {
  ok <- !is.na(y) & is.finite(p)
  y <- y[ok]; p <- p[ok]
  p <- pmin(pmax(p, 1e-8), 1 - 1e-8)
  lp <- stats::qlogis(p)
  i_fit <- try(stats::glm(y ~ 1 + offset(lp), family = stats::binomial()), silent = TRUE)
  s_fit <- try(stats::glm(y ~ lp, family = stats::binomial()), silent = TRUE)
  z <- stats::qnorm(1 - (1 - conf.level) / 2)
  if (inherits(i_fit, "try-error")) {
    intercept <- il <- ih <- NA_real_
  } else {
    intercept <- as.numeric(stats::coef(i_fit)[[1L]])
    se <- sqrt(diag(stats::vcov(i_fit)))[[1L]]
    il <- intercept - z * se; ih <- intercept + z * se
  }
  if (inherits(s_fit, "try-error") || length(stats::coef(s_fit)) < 2L) {
    slope <- sl <- sh <- NA_real_
  } else {
    slope <- as.numeric(stats::coef(s_fit)[[2L]])
    se <- sqrt(diag(stats::vcov(s_fit)))[[2L]]
    sl <- slope - z * se; sh <- slope + z * se
  }
  c(intercept = intercept, intercept_low = il, intercept_high = ih,
    slope = slope, slope_low = sl, slope_high = sh)
}

.r4vn_score_binary_performance <- function(y, p, seed = NULL, conf.level = 0.95) {
  auc <- .r4vn_score_auc_ci(y, p, conf.level = conf.level, seed = seed)
  cal <- .r4vn_score_calibration(y, p, conf.level = conf.level)
  ok <- !is.na(y) & is.finite(p)
  ee <- (y[ok] - p[ok])^2
  brier <- mean(ee)
  z <- stats::qnorm(1 - (1 - conf.level) / 2)
  bse <- if (length(ee) > 1L) stats::sd(ee) / sqrt(length(ee)) else NA_real_
  blo <- if (is.finite(bse)) max(0, brier - z * bse) else NA_real_
  bhi <- if (is.finite(bse)) min(1, brier + z * bse) else NA_real_
  c(AUC = auc[["est"]], AUC_low = auc[["low"]], AUC_high = auc[["high"]],
    Brier = brier, Brier_low = blo, Brier_high = bhi,
    Calibration_intercept = cal[["intercept"]],
    Calibration_intercept_low = cal[["intercept_low"]],
    Calibration_intercept_high = cal[["intercept_high"]],
    Calibration_slope = cal[["slope"]],
    Calibration_slope_low = cal[["slope_low"]],
    Calibration_slope_high = cal[["slope_high"]])
}

.r4vn_score_model_table <- function(fit, family, V = NULL) {
  b <- stats::coef(fit)
  if (is.null(V)) V <- stats::vcov(fit)
  V <- as.matrix(V)
  if (!all(names(b) %in% rownames(V)) || !all(names(b) %in% colnames(V))) V <- stats::vcov(fit)
  V <- V[names(b), names(b), drop = FALSE]
  se <- sqrt(diag(V))
  z <- b / se
  p <- 2 * stats::pnorm(abs(z), lower.tail = FALSE)
  lo <- b - stats::qnorm(.975) * se
  hi <- b + stats::qnorm(.975) * se
  keep <- names(b) != "(Intercept)"
  effect_name <- if (family == "logistic") "OR" else if (family == "cox") "HR" else "IRR/RR"
  data.frame(
    Term = names(b)[keep],
    Beta = unname(b[keep]),
    SE = unname(se[keep]),
    Effect = exp(unname(b[keep])),
    CI_low = exp(unname(lo[keep])),
    CI_high = exp(unname(hi[keep])),
    p = unname(p[keep]),
    Effect_measure = effect_name,
    stringsAsFactors = FALSE,
    check.names = FALSE
  )
}

.r4vn_score_model_publication <- function(fit, family, data, predictors, labels, V = NULL) {
  b <- stats::coef(fit)
  if (is.null(V)) V <- stats::vcov(fit)
  V <- as.matrix(V)
  if (!all(names(b) %in% rownames(V)) || !all(names(b) %in% colnames(V))) V <- stats::vcov(fit)
  V <- as.matrix(V)[names(b), names(b), drop = FALSE]
  tt <- stats::delete.response(stats::terms(fit))
  ref <- .r4vn_score_reference_row(data, predictors)
  xlev <- .r4vn_score_null(fit$xlevels, list())
  ctr <- fit$contrasts

  mmrow <- function(nd) {
    mf <- stats::model.frame(tt, data = nd, na.action = stats::na.pass, xlev = xlev)
    mm <- stats::model.matrix(tt, data = mf, contrasts.arg = ctr)
    z <- setNames(rep(0, length(b)), names(b))
    common <- intersect(colnames(mm), names(b))
    if (length(common)) z[common] <- as.numeric(mm[1L, common, drop = TRUE])
    z
  }
  xr <- mmrow(ref)
  effect_label <- if (family == "logistic") "OR (95% CI)" else if (family == "cox") "HR (95% CI)" else "IRR/RR (95% CI)"
  rows <- list(); k <- 0L
  for (v in predictors) {
    x <- data[[v]]
    lab <- labels[[v]]
    if (is.factor(x)) {
      levs <- levels(x)
      if (!length(levs)) next
      k <- k + 1L
      rows[[k]] <- data.frame(Predictor = lab, Category = levs[[1L]],
                              Effect = "1.00 (Reference)", p = "", stringsAsFactors = FALSE)
      if (length(levs) > 1L) {
        for (lev in levs[-1L]) {
          nd <- ref
          nd[[v]] <- factor(lev, levels = levs)
          cv <- mmrow(nd) - xr
          est <- sum(cv * b, na.rm = FALSE)
          vv <- as.numeric(t(cv) %*% V %*% cv)
          se <- if (is.finite(vv) && vv >= 0) sqrt(vv) else NA_real_
          lo <- est - stats::qnorm(.975) * se
          hi <- est + stats::qnorm(.975) * se
          pv <- if (is.finite(est) && is.finite(se) && se > 0) 2 * stats::pnorm(abs(est / se), lower.tail = FALSE) else NA_real_
          txt <- if (is.finite(est)) paste0(.r4vn_score_num(exp(est), 2), " (",
                                             .r4vn_score_num(exp(lo), 2), "\u2013",
                                             .r4vn_score_num(exp(hi), 2), ")") else ""
          ptxt <- if (!is.finite(pv)) "" else if (pv < .001) "<0.001" else .r4vn_score_num(pv, 3)
          k <- k + 1L
          rows[[k]] <- data.frame(Predictor = "", Category = lev, Effect = txt, p = ptxt,
                                  stringsAsFactors = FALSE)
        }
      }
    } else if (is.numeric(x)) {
      nd <- ref
      nd[[v]] <- as.numeric(nd[[v]]) + 1
      cv <- mmrow(nd) - xr
      est <- sum(cv * b, na.rm = FALSE)
      vv <- as.numeric(t(cv) %*% V %*% cv)
      se <- if (is.finite(vv) && vv >= 0) sqrt(vv) else NA_real_
      lo <- est - stats::qnorm(.975) * se
      hi <- est + stats::qnorm(.975) * se
      pv <- if (is.finite(est) && is.finite(se) && se > 0) 2 * stats::pnorm(abs(est / se), lower.tail = FALSE) else NA_real_
      txt <- if (is.finite(est)) paste0(.r4vn_score_num(exp(est), 2), " (",
                                         .r4vn_score_num(exp(lo), 2), "\u2013",
                                         .r4vn_score_num(exp(hi), 2), ")") else ""
      ptxt <- if (!is.finite(pv)) "" else if (pv < .001) "<0.001" else .r4vn_score_num(pv, 3)
      k <- k + 1L
      rows[[k]] <- data.frame(Predictor = lab, Category = "Per 1 unit", Effect = txt, p = ptxt,
                              stringsAsFactors = FALSE)
    }
  }
  if (!length(rows)) return(data.frame())
  out <- do.call(rbind, rows)
  names(out)[names(out) == "Effect"] <- effect_label
  rownames(out) <- NULL
  out
}

.r4vn_score_binom_ci <- function(x, n, conf.level = 0.95) {
  if (!is.finite(n) || n <= 0) return(c(est = NA, low = NA, high = NA))
  bt <- stats::binom.test(round(x), round(n), conf.level = conf.level)
  c(est = x / n, low = bt$conf.int[[1L]], high = bt$conf.int[[2L]])
}

.r4vn_score_cut_metrics <- function(y, positive) {
  ok <- !is.na(y) & !is.na(positive)
  y <- y[ok]; positive <- as.logical(positive[ok])
  TP <- sum(positive & y == 1); FN <- sum(!positive & y == 1)
  FP <- sum(positive & y == 0); TN <- sum(!positive & y == 0)
  sens <- .r4vn_score_binom_ci(TP, TP + FN)
  spec <- .r4vn_score_binom_ci(TN, TN + FP)
  ppv <- .r4vn_score_binom_ci(TP, TP + FP)
  npv <- .r4vn_score_binom_ci(TN, TN + FN)
  acc <- .r4vn_score_binom_ci(TP + TN, TP + TN + FP + FN)
  lrpos <- sens[["est"]] / pmax(1 - spec[["est"]], 1e-12)
  lrneg <- (1 - sens[["est"]]) / pmax(spec[["est"]], 1e-12)
  se_log_lrpos <- if (TP > 0 && TP + FN > 0 && FP > 0 && FP + TN > 0) sqrt(1 / TP - 1 / (TP + FN) + 1 / FP - 1 / (FP + TN)) else NA_real_
  se_log_lrneg <- if (FN > 0 && TP + FN > 0 && TN > 0 && FP + TN > 0) sqrt(1 / FN - 1 / (TP + FN) + 1 / TN - 1 / (FP + TN)) else NA_real_
  z <- stats::qnorm(.975)
  lrpos_ci <- if (is.finite(se_log_lrpos)) exp(log(lrpos) + c(-1, 1) * z * se_log_lrpos) else c(NA, NA)
  lrneg_ci <- if (is.finite(se_log_lrneg)) exp(log(lrneg) + c(-1, 1) * z * se_log_lrneg) else c(NA, NA)
  data.frame(
    Measure = c("Sensitivity", "Specificity", "PPV", "NPV", "Accuracy", "LR+", "LR\u2212"),
    Estimate = c(sens[1], spec[1], ppv[1], npv[1], acc[1], lrpos, lrneg),
    CI_low = c(sens[2], spec[2], ppv[2], npv[2], acc[2], lrpos_ci[1], lrneg_ci[1]),
    CI_high = c(sens[3], spec[3], ppv[3], npv[3], acc[3], lrpos_ci[2], lrneg_ci[2]),
    stringsAsFactors = FALSE
  )
}

.r4vn_score_roc_table <- function(y, score) {
  th <- sort(unique(score[is.finite(score)]))
  if (!length(th)) return(data.frame())
  do.call(rbind, lapply(th, function(t) {
    pos <- score >= t
    TP <- sum(pos & y == 1, na.rm = TRUE); FN <- sum(!pos & y == 1, na.rm = TRUE)
    FP <- sum(pos & y == 0, na.rm = TRUE); TN <- sum(!pos & y == 0, na.rm = TRUE)
    data.frame(score_cutoff = t,
               sensitivity = if ((TP + FN) > 0) TP / (TP + FN) else NA_real_,
               specificity = if ((TN + FP) > 0) TN / (TN + FP) else NA_real_,
               TP = TP, FN = FN, FP = FP, TN = TN)
  }))
}

.r4vn_score_find_risk_score <- function(risk_table, prob) {
  if (is.null(risk_table) || !nrow(risk_table) || is.null(prob) || !is.finite(prob)) return(NA_real_)
  rt <- risk_table[order(risk_table$Score), , drop = FALSE]
  idx <- which(rt$Predicted_risk >= prob)
  if (!length(idx)) NA_real_ else rt$Score[[idx[[1L]]]]
}

.r4vn_score_cutoff_table_binary <- function(y, score, risk_table = NULL, riskcut = NULL,
                                            sens_target = NULL, spec_target = NULL,
                                            cost_fp = 1, cost_fn = 1,
                                            cutoff_value = NULL, refprob = NULL, refcut = NULL) {
  roc <- .r4vn_score_roc_table(y, score)
  if (!nrow(roc)) return(data.frame())
  auc <- .r4vn_score_auc(y, score)
  roc$youden <- roc$sensitivity + roc$specificity - 1
  roc$iu <- abs(roc$sensitivity - auc) + abs(roc$specificity - auc)
  roc$cost <- cost_fn * roc$FN + cost_fp * roc$FP
  rows <- list()
  add <- function(method, r, target_prob = NA_real_) {
    data.frame(Method = method, Score_cutoff = r$score_cutoff, Risk_threshold = target_prob,
               Sensitivity = r$sensitivity, Specificity = r$specificity,
               Youden = r$youden, IU = r$iu, stringsAsFactors = FALSE)
  }
  rows[["youden"]] <- add("Youden", roc[which.max(roc$youden), , drop = FALSE])
  rows[["iu"]] <- add("IU", roc[which.min(roc$iu), , drop = FALSE])
  prev <- mean(y == 1, na.rm = TRUE)
  pc <- .r4vn_score_find_risk_score(risk_table, prev)
  if (is.finite(pc)) {
    rr <- roc[which.min(abs(roc$score_cutoff - pc)), , drop = FALSE]
    rows[["prevalence"]] <- add("Prevalence probability", rr, prev)
  }
  if (!is.null(riskcut) && length(riskcut)) {
    pc <- .r4vn_score_find_risk_score(risk_table, riskcut[[1L]])
    if (is.finite(pc)) {
      rr <- roc[which.min(abs(roc$score_cutoff - pc)), , drop = FALSE]
      rows[["risk"]] <- add("Clinical risk probability", rr, riskcut[[1L]])
    }
  }
  if (!is.null(sens_target) && is.finite(sens_target)) {
    z <- roc[roc$sensitivity >= sens_target, , drop = FALSE]
    if (nrow(z)) rows[["sens"]] <- add(paste0("Sensitivity \u2265 ", sens_target), z[which.max(z$specificity), , drop = FALSE])
  }
  if (!is.null(spec_target) && is.finite(spec_target)) {
    z <- roc[roc$specificity >= spec_target, , drop = FALSE]
    if (nrow(z)) rows[["spec"]] <- add(paste0("Specificity \u2265 ", spec_target), z[which.max(z$sensitivity), , drop = FALSE])
  }
  rows[["cost"]] <- add("Minimum misclassification cost", roc[which.min(roc$cost), , drop = FALSE])
  if (!is.null(cutoff_value) && is.finite(cutoff_value)) {
    rr <- roc[which.min(abs(roc$score_cutoff - cutoff_value)), , drop = FALSE]
    rows[["manual"]] <- add("Manual", rr)
  }
  if (!is.null(refprob) && !is.null(refcut) && is.finite(refcut)) {
    refclass <- as.integer(refprob >= refcut)
    r2 <- .r4vn_score_roc_table(refclass, score)
    if (nrow(r2)) {
      r2$youden <- r2$sensitivity + r2$specificity - 1
      rr0 <- r2[which.max(r2$youden), , drop = FALSE]
      # Report performance against the true outcome at the score threshold selected
      # to best reproduce the reference-probability classification.
      rr <- roc[which.min(abs(roc$score_cutoff - rr0$score_cutoff)), , drop = FALSE]
      rows[["refprob"]] <- add("Reference probability", rr, refcut)
    }
  }
  do.call(rbind, rows)
}

.r4vn_score_risk_group <- function(p, cuts) {
  if (is.null(cuts) || !length(cuts)) return(rep(NA_character_, length(p)))
  cuts <- sort(unique(cuts))
  k <- length(cuts) + 1L
  labs <- if (k == 2L) c("Low", "High") else if (k == 3L) c("Low", "Intermediate", "High") else if (k == 4L) c("Low", "Intermediate", "High", "Very high") else paste("Risk group", seq_len(k))
  as.character(cut(p, breaks = c(-Inf, cuts, Inf), labels = labs, right = FALSE))
}


.r4vn_score_add_observed_binary <- function(risk_table, y, score) {
  n <- table(score)
  ev <- tapply(y, score, sum, na.rm = TRUE)
  rr <- tapply(y, score, mean, na.rm = TRUE)
  key <- as.character(risk_table$Score)
  risk_table$Observed_n <- as.numeric(n[key])
  risk_table$Observed_events <- as.numeric(ev[key])
  risk_table$Observed_risk <- as.numeric(rr[key])
  risk_table$Observed_n[is.na(risk_table$Observed_n)] <- 0
  risk_table
}

.r4vn_score_add_observed_poisson <- function(risk_table, y, score) {
  n <- table(score)
  mn <- tapply(y, score, mean, na.rm = TRUE)
  key <- as.character(risk_table$Score)
  risk_table$Observed_n <- as.numeric(n[key])
  risk_table$Observed_mean <- as.numeric(mn[key])
  risk_table$Observed_n[is.na(risk_table$Observed_n)] <- 0
  risk_table
}

.r4vn_score_risk_table_binary <- function(y, score, range, riskcut = NULL, score_values = NULL) {
  fit <- .r4vn_score_recal_binary(y, score)
  seqs <- if (is.null(score_values)) seq(floor(range[["min"]]), ceiling(range[["max"]]), by = 1) else sort(unique(score_values))
  nd <- data.frame(.score = seqs)
  pr <- stats::predict(fit, newdata = nd, type = "link", se.fit = TRUE)
  z <- stats::qnorm(.975)
  p <- stats::plogis(pr$fit)
  lo <- stats::plogis(pr$fit - z * pr$se.fit)
  hi <- stats::plogis(pr$fit + z * pr$se.fit)
  out <- data.frame(Score = seqs, Predicted_risk = p, CI_low = lo, CI_high = hi,
                    Risk_group = .r4vn_score_risk_group(p, riskcut), stringsAsFactors = FALSE)
  list(table = out, fit = fit)
}

.r4vn_score_risk_table_poisson <- function(y, score, range, score_values = NULL) {
  fit <- stats::glm(y ~ score, family = stats::poisson(link = "log"))
  seqs <- if (is.null(score_values)) seq(floor(range[["min"]]), ceiling(range[["max"]]), by = 1) else sort(unique(score_values))
  nd <- data.frame(score = seqs)
  pr <- stats::predict(fit, newdata = nd, type = "link", se.fit = TRUE)
  z <- stats::qnorm(.975)
  out <- data.frame(Score = seqs, Predicted_mean = exp(pr$fit),
                    CI_low = exp(pr$fit - z * pr$se.fit), CI_high = exp(pr$fit + z * pr$se.fit),
                    stringsAsFactors = FALSE)
  list(table = out, fit = fit)
}

.r4vn_score_cox_risk <- function(fit, newdata, times) {
  bh <- survival::basehaz(fit, centered = FALSE)
  lp <- as.numeric(stats::predict(fit, newdata = newdata, type = "lp", reference = "zero"))
  H <- vapply(times, function(t) {
    j <- max(which(bh$time <= t), 0L)
    if (j == 0L) 0 else bh$hazard[[j]]
  }, numeric(1L))
  outer(exp(lp), H, function(e, h) 1 - exp(-h * e))
}

.r4vn_score_risk_table_cox <- function(time, event, score, range, times, riskcut = NULL, score_values = NULL) {
  if (is.null(times) || !length(times)) {
    times <- as.numeric(stats::quantile(time[event == 1], probs = c(.25, .5, .75), na.rm = TRUE, names = FALSE))
    times <- unique(times[is.finite(times) & times > 0])
  }
  d <- data.frame(.time = time, .event = event, .score = score)
  fit <- survival::coxph(survival::Surv(.time, .event) ~ .score, data = d, x = TRUE, y = TRUE)
  seqs <- if (is.null(score_values)) seq(floor(range[["min"]]), ceiling(range[["max"]]), by = 1) else sort(unique(score_values))
  out <- list()
  for (s in seqs) {
    sf <- try(survival::survfit(fit, newdata = data.frame(.score = s)), silent = TRUE)
    if (inherits(sf, "try-error")) next
    su <- summary(sf, times = times, extend = TRUE)
    surv <- as.numeric(su$surv); low_s <- as.numeric(su$lower); high_s <- as.numeric(su$upper)
    out[[length(out) + 1L]] <- data.frame(
      Score = s, Time = times,
      Predicted_risk = 1 - surv,
      CI_low = 1 - high_s,
      CI_high = 1 - low_s,
      stringsAsFactors = FALSE
    )
  }
  tab <- do.call(rbind, out)
  if (!is.null(riskcut) && length(riskcut)) tab$Risk_group <- .r4vn_score_risk_group(tab$Predicted_risk, riskcut)
  list(table = tab, fit = fit, times = times)
}

.r4vn_score_km_censor <- function(time, event) {
  survival::survfit(survival::Surv(time, 1 - event) ~ 1)
}

.r4vn_score_km_value <- function(sf, t, left = FALSE) {
  if (!length(sf$time)) return(rep(1, length(t)))
  eps <- sqrt(.Machine$double.eps) * pmax(1, abs(t))
  tt <- if (left) t - eps else t
  idx <- findInterval(tt, sf$time)
  out <- rep(1, length(tt))
  pos <- idx > 0L
  out[pos] <- sf$surv[idx[pos]]
  pmax(out, 1e-6)
}

.r4vn_score_time_roc <- function(time, event, marker, horizon) {
  sfG <- .r4vn_score_km_censor(time, event)
  case <- event == 1 & time <= horizon
  control <- time > horizon
  wc <- rep(0, length(time)); wn <- rep(0, length(time))
  wc[case] <- 1 / .r4vn_score_km_value(sfG, time[case], left = TRUE)
  wn[control] <- 1 / .r4vn_score_km_value(sfG, rep(horizon, sum(control)), left = FALSE)
  th <- sort(unique(marker[is.finite(marker)]))
  roc <- do.call(rbind, lapply(th, function(cut) {
    pos <- marker >= cut
    sens <- sum(wc * pos, na.rm = TRUE) / pmax(sum(wc), 1e-12)
    spec <- sum(wn * (!pos), na.rm = TRUE) / pmax(sum(wn), 1e-12)
    TP <- sum(wc * pos, na.rm = TRUE); FN <- sum(wc * (!pos), na.rm = TRUE)
    FP <- sum(wn * pos, na.rm = TRUE); TN <- sum(wn * (!pos), na.rm = TRUE)
    data.frame(score_cutoff = cut, sensitivity = sens, specificity = spec, TP = TP, FN = FN, FP = FP, TN = TN)
  }))
  # Add ROC endpoints and integrate the cumulative/dynamic IPCW ROC curve.
  xy <- rbind(data.frame(fpr = 0, tpr = 0),
              data.frame(fpr = 1 - roc$specificity, tpr = roc$sensitivity),
              data.frame(fpr = 1, tpr = 1))
  xy <- xy[order(xy$fpr, xy$tpr), , drop = FALSE]
  auc <- sum(diff(xy$fpr) * (head(xy$tpr, -1L) + tail(xy$tpr, -1L)) / 2)
  roc$youden <- roc$sensitivity + roc$specificity - 1
  roc$iu <- abs(roc$sensitivity - auc) + abs(roc$specificity - auc)
  attr(roc, "auc") <- auc
  attr(roc, "weights") <- list(case = wc, control = wn)
  roc
}

.r4vn_score_cutoff_table_cox <- function(time, event, score, horizon, risk_table = NULL,
                                         riskcut = NULL, sens_target = NULL, spec_target = NULL,
                                         cost_fp = 1, cost_fn = 1, cutoff_value = NULL) {
  roc <- .r4vn_score_time_roc(time, event, score, horizon)
  auc <- attr(roc, "auc")
  roc$cost <- cost_fn * roc$FN + cost_fp * roc$FP
  rows <- list()
  add <- function(method, r, target_prob = NA_real_) data.frame(
    Method = method, Score_cutoff = r$score_cutoff, Risk_threshold = target_prob,
    Sensitivity = r$sensitivity, Specificity = r$specificity,
    Youden = r$youden, IU = r$iu, Time = horizon, AUC = auc, stringsAsFactors = FALSE)
  rows[["youden"]] <- add("Time-dependent Youden (IPCW)", roc[which.max(roc$youden), , drop = FALSE])
  rows[["iu"]] <- add("Time-dependent IU (IPCW)", roc[which.min(roc$iu), , drop = FALSE])
  if (!is.null(riskcut) && length(riskcut) && !is.null(risk_table)) {
    rt <- risk_table[risk_table$Time == horizon, , drop = FALSE]
    pc <- .r4vn_score_find_risk_score(rt, riskcut[[1L]])
    if (is.finite(pc)) {
      rr <- roc[which.min(abs(roc$score_cutoff - pc)), , drop = FALSE]
      rows[["risk"]] <- add("Clinical risk probability", rr, riskcut[[1L]])
    }
  }
  if (!is.null(sens_target) && is.finite(sens_target)) {
    z <- roc[roc$sensitivity >= sens_target, , drop = FALSE]
    if (nrow(z)) rows[["sens"]] <- add(paste0("Sensitivity \u2265 ", sens_target), z[which.max(z$specificity), , drop = FALSE])
  }
  if (!is.null(spec_target) && is.finite(spec_target)) {
    z <- roc[roc$specificity >= spec_target, , drop = FALSE]
    if (nrow(z)) rows[["spec"]] <- add(paste0("Specificity \u2265 ", spec_target), z[which.max(z$sensitivity), , drop = FALSE])
  }
  rows[["cost"]] <- add("Minimum IPCW misclassification cost", roc[which.min(roc$cost), , drop = FALSE])
  if (!is.null(cutoff_value) && is.finite(cutoff_value)) {
    rr <- roc[which.min(abs(roc$score_cutoff - cutoff_value)), , drop = FALSE]
    rows[["manual"]] <- add("Manual", rr)
  }
  do.call(rbind, rows)
}

.r4vn_score_cut_metrics_cox <- function(time, event, score, cutoff, horizon) {
  roc <- .r4vn_score_time_roc(time, event, score, horizon)
  r <- roc[which.min(abs(roc$score_cutoff - cutoff)), , drop = FALSE]
  TP <- r$TP; FN <- r$FN; FP <- r$FP; TN <- r$TN
  sens <- TP / pmax(TP + FN, 1e-12); spec <- TN / pmax(TN + FP, 1e-12)
  ppv <- TP / pmax(TP + FP, 1e-12); npv <- TN / pmax(TN + FN, 1e-12)
  acc <- (TP + TN) / pmax(TP + TN + FP + FN, 1e-12)
  lrpos <- sens / pmax(1 - spec, 1e-12); lrneg <- (1 - sens) / pmax(spec, 1e-12)
  data.frame(Measure = c("Sensitivity (IPCW)", "Specificity (IPCW)", "PPV (IPCW)", "NPV (IPCW)",
                         "Accuracy (IPCW)", "LR+ (IPCW)", "LR\u2212 (IPCW)"),
             Estimate = c(sens, spec, ppv, npv, acc, lrpos, lrneg),
             CI_low = NA_real_, CI_high = NA_real_, Time = horizon, stringsAsFactors = FALSE)
}

.r4vn_score_ipcw_brier <- function(time, event, pred_risk, horizon) {
  sfG <- .r4vn_score_km_censor(time, event)
  case <- event == 1 & time <= horizon
  control <- time > horizon
  w <- rep(0, length(time))
  w[case] <- 1 / .r4vn_score_km_value(sfG, time[case], left = TRUE)
  w[control] <- 1 / .r4vn_score_km_value(sfG, rep(horizon, sum(control)), left = FALSE)
  surv_obs <- as.numeric(time > horizon)
  pred_surv <- 1 - pred_risk
  sum(w * (surv_obs - pred_surv)^2, na.rm = TRUE) / length(time)
}

.r4vn_score_reference_probability <- function(refprob, data, env) {
  if (is.null(refprob)) return(NULL)
  if (is.character(refprob) && length(refprob) == 1L && refprob %in% names(data)) return(as.numeric(data[[refprob]]))
  if (is.symbol(substitute(refprob))) {
    nm <- as.character(substitute(refprob))
    if (nm %in% names(data)) return(as.numeric(data[[nm]]))
  }
  val <- try(eval(substitute(refprob), env), silent = TRUE)
  if (!inherits(val, "try-error") && is.numeric(val) && length(val) == nrow(data)) return(as.numeric(val))
  if (is.numeric(refprob) && length(refprob) == nrow(data)) return(as.numeric(refprob))
  stop("refprob must be a numeric vector with one value per row or the name of a probability variable in data.", call. = FALSE)
}

.r4vn_score_refprob_table <- function(refprob, scoreprob) {
  if (is.null(refprob)) return(NULL)
  ok <- is.finite(refprob) & is.finite(scoreprob)
  d <- refprob[ok] - scoreprob[ok]
  data.frame(
    Measure = c("Correlation", "MAE", "RMSE", "Mean difference (reference \u2212 score)"),
    Estimate = c(stats::cor(refprob[ok], scoreprob[ok]), mean(abs(d)), sqrt(mean(d^2)), mean(d)),
    stringsAsFactors = FALSE
  )
}

.r4vn_score_auc_compare <- function(y, p1, p2, seed = NULL, B = 500L) {
  a1 <- .r4vn_score_auc(y, p1); a2 <- .r4vn_score_auc(y, p2); delta <- a2 - a1
  if (requireNamespace("pROC", quietly = TRUE)) {
    r1 <- try(pROC::roc(y, p1, quiet = TRUE, direction = "<"), silent = TRUE)
    r2 <- try(pROC::roc(y, p2, quiet = TRUE, direction = "<"), silent = TRUE)
    if (!inherits(r1, "try-error") && !inherits(r2, "try-error")) {
      tt <- try(pROC::roc.test(r1, r2, paired = TRUE, method = "delong"), silent = TRUE)
      if (!inherits(tt, "try-error")) {
        ci <- try(as.numeric(tt$conf.int), silent = TRUE)
        okci <- !inherits(ci, "try-error") && length(ci) >= 2L
        return(c(delta = delta, low = if (okci) ci[[1L]] else NA,
                 high = if (okci) ci[[2L]] else NA, p = as.numeric(tt$p.value)))
      }
    }
  }
  .r4vn_set_seed_if(seed)
  n <- length(y)
  dd <- replicate(B, {
    ii <- sample.int(n, n, replace = TRUE)
    .r4vn_score_auc(y[ii], p2[ii]) - .r4vn_score_auc(y[ii], p1[ii])
  })
  dd <- dd[is.finite(dd)]
  if (!length(dd)) return(c(delta = delta, low = NA_real_, high = NA_real_, p = NA_real_))
  ci <- stats::quantile(dd, c(.025, .975), na.rm = TRUE, names = FALSE)
  p <- min(1, 2 * min(mean(dd <= 0), mean(dd >= 0)))
  c(delta = delta, low = ci[[1L]], high = ci[[2L]], p = p)
}

.r4vn_score_comparison_binary <- function(y, pred_original, pred_score_model, pred_clinical, seed = NULL) {
  vals <- list(`Original model` = pred_original, `Model score` = pred_score_model, `Clinical score` = pred_clinical)
  tab <- do.call(rbind, lapply(names(vals), function(nm) {
    m <- .r4vn_score_binary_performance(y, vals[[nm]], seed = seed)
    data.frame(Model = nm, AUC = m[["AUC"]], AUC_low = m[["AUC_low"]], AUC_high = m[["AUC_high"]],
               Brier = m[["Brier"]], Brier_low = m[["Brier_low"]], Brier_high = m[["Brier_high"]],
               Calibration_intercept = m[["Calibration_intercept"]],
               Calibration_intercept_low = m[["Calibration_intercept_low"]],
               Calibration_intercept_high = m[["Calibration_intercept_high"]],
               Calibration_slope = m[["Calibration_slope"]],
               Calibration_slope_low = m[["Calibration_slope_low"]],
               Calibration_slope_high = m[["Calibration_slope_high"]], stringsAsFactors = FALSE)
  }))
  base_auc <- tab$AUC[tab$Model == "Original model"]
  base_brier <- tab$Brier[tab$Model == "Original model"]
  tab$Delta_AUC_vs_original <- tab$AUC - base_auc
  tab$Delta_AUC_low <- NA_real_
  tab$Delta_AUC_high <- NA_real_
  tab$Delta_AUC_p <- NA_real_
  tab$Delta_Brier_vs_original <- tab$Brier - base_brier
  for (nm in c("Model score", "Clinical score")) {
    j <- which(tab$Model == nm)
    pp <- if (nm == "Model score") pred_score_model else pred_clinical
    cmp <- .r4vn_score_auc_compare(y, pred_original, pp, seed = seed)
    tab$Delta_AUC_low[j] <- cmp[["low"]]
    tab$Delta_AUC_high[j] <- cmp[["high"]]
    tab$Delta_AUC_p[j] <- cmp[["p"]]
  }
  tab
}

.r4vn_score_comparison_cox <- function(time, event, original_fit, score_model_fit, score_fit,
                                       raw_data, score_data, score, times) {
  lp0 <- as.numeric(stats::predict(original_fit, newdata = raw_data, type = "lp"))
  lp1 <- as.numeric(stats::predict(score_model_fit, newdata = score_data, type = "lp"))
  lp2 <- score
  c0 <- .r4vn_score_cindex(time, event, lp0); c1 <- .r4vn_score_cindex(time, event, lp1); c2 <- .r4vn_score_cindex(time, event, lp2)
  tab <- data.frame(Model = c("Original model", "Model score", "Clinical score"),
                    C_index = c(c0[1], c1[1], c2[1]), C_low = c(c0[2], c1[2], c2[2]), C_high = c(c0[3], c1[3], c2[3]),
                    stringsAsFactors = FALSE)
  tab$Delta_C_index_vs_original <- tab$C_index - tab$C_index[[1L]]
  if (!is.null(times) && length(times)) {
    pr0 <- .r4vn_score_cox_risk(original_fit, raw_data, times)
    pr1 <- .r4vn_score_cox_risk(score_model_fit, score_data, times)
    d2 <- data.frame(.r4vn_score_time = time, .r4vn_score_y = event, .score = score)
    sf2 <- survival::coxph(survival::Surv(.r4vn_score_time, .r4vn_score_y) ~ .score, data = d2, x = TRUE, y = TRUE)
    pr2 <- .r4vn_score_cox_risk(sf2, data.frame(.score = score), times)
    b <- do.call(rbind, lapply(seq_along(times), function(j) data.frame(
      Time = times[[j]],
      Original_model = .r4vn_score_ipcw_brier(time, event, pr0[, j], times[[j]]),
      Model_score = .r4vn_score_ipcw_brier(time, event, pr1[, j], times[[j]]),
      Clinical_score = .r4vn_score_ipcw_brier(time, event, pr2[, j], times[[j]])
    )))
    attr(tab, "time_brier") <- b
  }
  tab
}

.r4vn_score_poisson_perf <- function(y, p) {
  data.frame(RMSE = sqrt(mean((y - p)^2)), MAE = mean(abs(y - p)), Mean_prediction = mean(p), stringsAsFactors = FALSE)
}

.r4vn_score_decision_curve <- function(y, pred_original, pred_clinical, thresholds = seq(.01, .50, by = .01)) {
  n <- length(y); prev <- mean(y == 1)
  out <- lapply(thresholds, function(pt) {
    nb <- function(p) {
      pos <- p >= pt
      TP <- sum(pos & y == 1); FP <- sum(pos & y == 0)
      TP / n - FP / n * pt / (1 - pt)
    }
    data.frame(Threshold = pt, Original_model = nb(pred_original), Clinical_score = nb(pred_clinical),
               Treat_all = prev - (1 - prev) * pt / (1 - pt), Treat_none = 0)
  })
  do.call(rbind, out)
}

.r4vn_score_calibration_data <- function(y, p, groups = 10L) {
  br <- unique(stats::quantile(p, probs = seq(0, 1, length.out = groups + 1L), na.rm = TRUE, names = FALSE))
  if (length(br) < 3L) return(data.frame(Predicted = mean(p), Observed = mean(y), n = length(y)))
  g <- cut(p, breaks = br, include.lowest = TRUE)
  do.call(rbind, lapply(split(seq_along(y), g), function(ii) data.frame(Predicted = mean(p[ii]), Observed = mean(y[ii]), n = length(ii))))
}

.r4vn_score_orientation_display <- function(dict, effect_measure = "Effect ratio") {
  if (is.null(dict) || !nrow(dict)) return(data.frame())
  original_name <- paste0("Model-reference ", effect_measure)
  risk_name <- paste0("Risk-oriented ", effect_measure)
  out <- data.frame(
    Predictor = dict$predictor_label,
    Category = dict$category,
    original = ifelse(dict$model_reference, "1.00 (Model ref.)", .r4vn_score_num(dict$model_ratio, 2)),
    risk = ifelse(dict$scoring_reference, "1.00 (Score ref.)", .r4vn_score_num(dict$scoring_ratio, 2)),
    Point = as.character(dict$clinical_point),
    stringsAsFactors = FALSE,
    check.names = FALSE
  )
  names(out)[3:4] <- c(original_name, risk_name)
  dup <- duplicated(out$Predictor)
  out$Predictor[dup] <- ""
  rownames(out) <- NULL
  out
}

.r4vn_score_scorecard_display <- function(dict, range) {
  out <- data.frame(Predictor = dict$predictor_label, Category = dict$category,
                    Point = as.character(dict$clinical_point), stringsAsFactors = FALSE)
  dup <- duplicated(out$Predictor)
  out$Predictor[dup] <- ""
  out <- rbind(out, data.frame(Predictor = "Total score", Category = "",
                               Point = paste0(range[["min"]], "\u2013", range[["max"]]), stringsAsFactors = FALSE))
  rownames(out) <- NULL
  out
}

.r4vn_score_risk_display <- function(tab, family) {
  if (is.null(tab) || !nrow(tab)) return(tab)
  out <- tab
  if (family %in% c("logistic", "cox") || "Predicted_risk" %in% names(out)) {
    for (nm in intersect(c("Predicted_risk", "CI_low", "CI_high"), names(out))) out[[nm]] <- .r4vn_score_pct(out[[nm]], 1)
  } else {
    for (nm in intersect(c("Predicted_mean", "CI_low", "CI_high"), names(out))) out[[nm]] <- .r4vn_score_num(out[[nm]], 3)
  }
  names(out) <- gsub("_", " ", names(out), fixed = TRUE)
  out
}

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

.r4vn_score_html_table <- function(df, title = NULL, scorecard = FALSE) {
  if (is.null(df) || !nrow(df)) return("")
  d <- df
  d[] <- lapply(d, function(z) {
    if (is.numeric(z)) {
      ifelse(is.na(z), "", formatC(z, format = "f", digits = 3))
    } else ifelse(is.na(z), "", as.character(z))
  })
  head <- paste0("<th>", .r4vn_score_html_escape(names(d)), "</th>", collapse = "")
  body <- character(nrow(d))
  for (i in seq_len(nrow(d))) {
    cls <- if (scorecard && identical(as.character(d[i, 1L]), "Total score")) " class='total'" else ""
    cells <- paste0("<td>", .r4vn_score_html_escape(unlist(d[i, , drop = FALSE], use.names = FALSE)), "</td>", collapse = "")
    body[[i]] <- paste0("<tr", cls, ">", cells, "</tr>")
  }
  paste0(
    if (!is.null(title)) paste0("<h3>", .r4vn_score_html_escape(title), "</h3>") else "",
    "<div class='table-wrap'><table", if (scorecard) " class='scorecard'" else "", "><thead><tr>", head, "</tr></thead><tbody>",
    paste(body, collapse = ""), "</tbody></table></div>"
  )
}


.r4vn_score_view_transpose <- function(df, id = "Model", max_rows = 4L, min_cols = 9L) {
  if (is.null(df) || !is.data.frame(df) || !nrow(df) || !id %in% names(df) ||
      nrow(df) > max_rows || ncol(df) < min_cols) return(df)
  labs <- as.character(df[[id]])
  metrics <- setdiff(names(df), id)
  out <- data.frame(Measure = gsub("_", " ", metrics, fixed = TRUE),
                    stringsAsFactors = FALSE, check.names = FALSE)
  format_one <- function(z, nm) {
    if (is.numeric(z)) {
      if (is.na(z)) return("")
      if (grepl("(^|_)p($|_)|p.value|p value", nm, ignore.case = TRUE)) {
        if (z < .001) return("<0.001")
        return(.r4vn_score_num(z, 3))
      }
      return(.r4vn_score_num(z, 3))
    }
    if (is.na(z)) "" else as.character(z)
  }
  for (i in seq_along(labs)) {
    nm <- labs[[i]]
    if (is.na(nm) || !nzchar(nm)) nm <- paste0("Model ", i)
    out[[nm]] <- vapply(metrics, function(m) format_one(df[[m]][[i]], m), character(1))
  }
  out
}

.r4vn_score_available_plots <- function(x) {
  pd <- x$plot_data
  if (is.null(pd) || !length(pd)) return(character())
  ord <- c("risk", "roc", "calibration", "decision", "distribution")
  ord[vapply(ord, function(nm) {
    z <- pd[[nm]]
    is.data.frame(z) && nrow(z) > 0L
  }, logical(1))]
}

.r4vn_score_plot_title <- function(x, which) {
  switch(
    which,
    risk = if (identical(x$family, "poisson") && !isTRUE(x$binary_outcome))
      "Score-to-expected-value curve" else "Score-to-risk curve",
    roc = if (identical(x$family, "cox") && !is.null(x$cutoff_time))
      paste0("Time-dependent ROC curve at time ", x$cutoff_time) else
      "Receiver operating characteristic curve",
    calibration = "Calibration plot",
    decision = "Decision-curve analysis",
    distribution = if (identical(x$family, "poisson") && !isTRUE(x$binary_outcome))
      "Clinical score and observed count" else "Clinical score distribution by outcome",
    tools::toTitleCase(which)
  )
}

.r4vn_score_plot_draw <- function(x, which, title = NULL, font_family = "sans") {
  d <- x$plot_data[[which]]
  if (is.null(d) || !is.data.frame(d) || !nrow(d))
    stop("No plot data are available for '", which, "'.", call. = FALSE)
  old <- graphics::par(no.readonly = TRUE)
  on.exit(graphics::par(old), add = TRUE)
  graphics::par(family = font_family)
  main <- .r4vn_score_null(title, .r4vn_score_plot_title(x, which))

  if (which == "risk") {
    yname <- if ("Predicted_risk" %in% names(d)) "Predicted_risk" else "Predicted_mean"
    if ("Time" %in% names(d)) {
      tt <- unique(d$Time)
      yr <- range(d[[yname]], na.rm = TRUE)
      xr <- range(d$Score, na.rm = TRUE)
      if (diff(xr) == 0) xr <- xr + c(-.5, .5)
      if (diff(yr) == 0) yr <- yr + c(-.05, .05) * max(1, abs(yr[[1L]]))
      graphics::plot(NA, xlim = xr, ylim = yr, xlab = "Score",
                     ylab = if (yname == "Predicted_risk") "Predicted risk" else "Predicted value",
                     main = main)
      for (i in seq_along(tt)) {
        z <- d[d$Time == tt[[i]], , drop = FALSE]
        graphics::lines(z$Score, z[[yname]], type = "b", lty = i, pch = i)
      }
      graphics::legend("topleft", legend = paste0("Time ", tt),
                       lty = seq_along(tt), pch = seq_along(tt), bty = "n")
    } else {
      graphics::plot(d$Score, d[[yname]], type = "b", xlab = "Score",
                     ylab = if (yname == "Predicted_risk") "Predicted risk" else "Predicted value",
                     main = main)
    }
  } else if (which == "roc") {
    fpr <- pmax(0, pmin(1, 1 - d$specificity))
    tpr <- pmax(0, pmin(1, d$sensitivity))
    graphics::plot(fpr, tpr, type = "l", lwd = 2,
                   xlim = c(0, 1), ylim = c(0, 1), xaxs = "i", yaxs = "i",
                   xlab = "1 - Specificity", ylab = "Sensitivity", asp = 1,
                   main = main)
    graphics::abline(0, 1, lty = 2, col = "gray60")
    if (!is.null(x$selected_cutoff) && is.finite(x$selected_cutoff) &&
        "score_cutoff" %in% names(d)) {
      j <- which.min(abs(d$score_cutoff - x$selected_cutoff))
      graphics::points(fpr[[j]], tpr[[j]], pch = 19)
      graphics::text(fpr[[j]], tpr[[j]],
                     labels = paste0("cutoff >=", x$selected_cutoff), pos = 4, cex = .85)
    }
  } else if (which == "calibration") {
    graphics::plot(d$Predicted, d$Observed, type = "b",
                   xlim = c(0, 1), ylim = c(0, 1), xaxs = "i", yaxs = "i",
                   xlab = "Predicted risk", ylab = "Observed risk", asp = 1,
                   main = main)
    graphics::abline(0, 1, lty = 2, col = "gray60")
  } else if (which == "decision") {
    yy <- cbind(d$Original_model, d$Clinical_score, d$Treat_all, d$Treat_none)
    graphics::matplot(d$Threshold, yy, type = "l", lty = 1:4, lwd = c(2, 2, 1, 1),
                      xlab = "Threshold probability", ylab = "Net benefit",
                      main = main)
    graphics::abline(h = 0, col = "gray80")
    graphics::legend("topright",
                     legend = c("Original model", "Clinical score", "Treat all", "Treat none"),
                     lty = 1:4, lwd = c(2, 2, 1, 1), bty = "n", cex = .85)
  } else if (which == "distribution") {
    if (identical(x$family, "poisson") && !isTRUE(x$binary_outcome)) {
      graphics::plot(d$Score, d$Outcome,
                     xlab = "Clinical score", ylab = "Observed count", main = main)
    } else {
      dd <- data.frame(
        Score = d$Score,
        Outcome = factor(d$Outcome, levels = c(0, 1), labels = c("No event", "Event"))
      )
      graphics::boxplot(Score ~ Outcome, data = dd,
                        xlab = "Outcome", ylab = "Clinical score", main = main)
    }
  }
  invisible(d)
}

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

.r4vn_score_plot_png <- function(x, which, width = 1152L, height = 792L, res = 144L) {
  path <- tempfile(fileext = ".png")
  on.exit(unlink(path), add = TRUE)
  args <- list(filename = path, width = width, height = height, units = "px",
               res = res, bg = "white")
  if (isTRUE(capabilities("cairo"))) args$type <- "cairo-png"
  do.call(grDevices::png, args)
  tryCatch(.r4vn_score_plot_draw(x, which), finally = grDevices::dev.off())
  size <- file.info(path)$size
  if (!is.finite(size) || size <= 0) stop("Could not render plot as PNG.", call. = FALSE)
  bytes <- readBin(path, what = "raw", n = size)
  paste0(
    "<img class='score-plot-image' alt='",
    .r4vn_score_html_escape(.r4vn_score_plot_title(x, which)),
    "' src='data:image/png;base64,", .r4vn_score_base64(bytes), "'>"
  )
}

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

.r4vn_score_plot_html <- function(x) {
  if (!isTRUE(x$plot_enabled)) return("")
  available <- .r4vn_score_available_plots(x)
  if (!length(available)) return("")
  blocks <- vapply(seq_along(available), function(i) {
    nm <- available[[i]]
    fig <- try(.r4vn_score_plot_png(x, nm), silent = TRUE)
    if (inherits(fig, "try-error")) fig <- try(.r4vn_score_plot_svg(x, nm), silent = TRUE)
    if (inherits(fig, "try-error")) return("")
    paste0(
      "<section class='figure'><h2>Figure ", i, ". ",
      .r4vn_score_html_escape(.r4vn_score_plot_title(x, nm)),
      "</h2><div class='score-chart'>", fig, "</div></section>"
    )
  }, character(1))
  paste0(blocks[nzchar(blocks)], collapse = "")
}

.r4vn_score_show_html <- function(x, show = TRUE) {
  path <- tempfile("r4vn-tabscore-", fileext = ".html")
  comparison_view <- .r4vn_score_view_transpose(x$tables$comparison)
  external_view <- .r4vn_score_view_transpose(x$tables$external_validation)
  plot_html <- .r4vn_score_plot_html(x)
  css <- "
  body{font-family:Arial,Helvetica,sans-serif;margin:26px;color:#111;background:#fff;line-height:1.35}
  h1{font-size:24px;margin:0 0 4px 0} h2{font-size:19px;margin-top:30px;border-bottom:1px solid #ddd;padding-bottom:6px}
  h3{font-size:16px;margin:22px 0 8px 0}.meta{color:#555;margin-bottom:18px}.note{background:#f6f6f6;border-left:4px solid #aaa;padding:10px 12px;margin:14px 0}
  .table-wrap{overflow-x:auto;margin-bottom:18px}table{border-collapse:collapse;width:100%;font-size:14px}th{text-align:left;border-bottom:1px solid #aaa;padding:8px 10px;white-space:nowrap}
  td{border-bottom:1px solid #e6e6e6;padding:8px 10px;vertical-align:top}tr.total td{font-weight:700;border-top:1.5px solid #777;border-bottom:1.5px solid #777}
  table.scorecard th:last-child,table.scorecard td:last-child{text-align:right}table.scorecard th:first-child{width:34%}table.scorecard th:nth-child(2){width:46%}
  .small{font-size:12px;color:#666}.section{margin-bottom:24px}.figure{margin:30px 0 34px 0}.score-chart{max-width:920px;overflow:hidden;margin:12px 0 0 0}
  .score-chart img.score-plot-image,.score-chart svg{display:block;width:100%;height:auto;max-width:920px;background:#fff}.score-chart svg text{font-family:Arial,Helvetica,sans-serif!important;letter-spacing:normal!important;word-spacing:normal!important}"
  parts <- c(
    "<!doctype html><html><head><meta charset='utf-8'><meta name='viewport' content='width=device-width,initial-scale=1'><title>R4VN Scorecard</title><style>", css, "</style></head><body>",
    "<h1>R4VN Scorecard</h1>",
    paste0("<div class='meta'>Model: <b>", .r4vn_score_html_escape(x$family),
           "</b> &nbsp; | &nbsp; n = ", x$n,
           " &nbsp; | &nbsp; Predictors: ", length(x$selected_predictors), "</div>"),
    paste0("<div class='note'><b>Clinical score:</b> theoretical range ", x$score_range[["min"]], "\u2013", x$score_range[["max"]],
           "; observed development range ", x$observed_score_range[[1L]], "\u2013", x$observed_score_range[[2L]],
           if (!is.null(x$selected_cutoff)) paste0("; selected cutoff \u2265", x$selected_cutoff, " (", .r4vn_score_html_escape(x$selected_cutoff_method), ")") else "",
           ".</div>"),
    if (!is.null(x$representation_warning)) paste0("<div class='note'><b>Predictor simplification note:</b> ", .r4vn_score_html_escape(x$representation_warning), "</div>") else "",
    if (!is.null(x$simplification_warning)) paste0("<div class='note'><b>Point simplification note:</b> ", .r4vn_score_html_escape(x$simplification_warning), "</div>") else "",
    .r4vn_score_html_table(x$tables$model, "1. Final prediction model"),
    .r4vn_score_html_table(x$tables$risk_orientation, "2. Risk-oriented coding for the score"),
    .r4vn_score_html_table(x$tables$scorecard, "3. Clinical scorecard", scorecard = TRUE),
    paste0("<p class='small'>Total score is the sum of one non-negative risk contribution per predictor when riskonly=TRUE. A protective contrast in the fitted model is automatically re-referenced to the lowest-risk category for scoring; the original fitted-model reference and estimates are preserved in the final-model table.</p>"),
    .r4vn_score_html_table(x$tables$risk, "4. Score-to-risk conversion"),
    .r4vn_score_html_table(x$tables$cutoff, "5. Cutoff selection"),
    .r4vn_score_html_table(x$tables$cutoff_performance, "6. Performance at the selected cutoff"),
    .r4vn_score_html_table(comparison_view, "7. Original model versus scorecard"),
    .r4vn_score_html_table(x$tables$time_brier, "Time-specific IPCW Brier scores"),
    .r4vn_score_html_table(x$tables$reference_probability, "Reference probability comparison"),
    .r4vn_score_html_table(x$tables$validation, "Internal validation"),
    .r4vn_score_html_table(external_view, "External validation"),
    plot_html,
    "</body></html>"
  )
  writeLines(enc2utf8(parts), path, useBytes = TRUE)
  if (isTRUE(show)) {
    viewer <- getOption("viewer")
    if (is.function(viewer)) viewer(path) else utils::browseURL(path)
  }
  normalizePath(path, winslash = "/", mustWork = TRUE)
}

.r4vn_score_bootstrap_validation <- function(object, B = 500L, seed = NULL) {
  if (B < 20L) warning("bootstrap < 20 is useful only for quick testing, not final validation.", call. = FALSE)
  args <- object$.rebuild_args
  dat <- object$.development_data_raw
  n <- nrow(dat)
  .r4vn_set_seed_if(seed)
  out <- vector("list", B)
  for (b in seq_len(B)) {
    ii <- sample.int(n, n, replace = TRUE)
    db <- dat[ii, , drop = FALSE]
    a <- args
    a$data <- db
    a$validate <- "none"
    a$show <- FALSE
    a$console <- FALSE
    a$plot <- FALSE
    # Avoid nested AUC-CI bootstraps while validating the whole pipeline.
    # Predictions are always retained, so comparison tables are unnecessary here.
    a$compare <- FALSE
    a$decision <- FALSE
    a$bootstrap <- 0L
    fitb <- try(do.call(tabscore, a), silent = TRUE)
    if (inherits(fitb, "try-error")) next
    if (object$family == "logistic" || (object$family == "poisson" && object$binary_outcome)) {
      yb <- fitb$.development_data$.r4vn_score_y
      pb_app_o <- fitb$predictions$original
      pb_app_s <- fitb$predictions$clinical
      po_test <- try(stats::predict(fitb, newdata = dat, type = "model"), silent = TRUE)
      ps_test <- try(stats::predict(fitb, newdata = dat, type = "risk"), silent = TRUE)
      score_test <- try(stats::predict(fitb, newdata = dat, type = "score"), silent = TRUE)
      if (inherits(po_test, "try-error") || inherits(ps_test, "try-error") || inherits(score_test, "try-error")) next
      yt <- object$.development_data$.r4vn_score_y
      co_a <- .r4vn_score_calibration(yb, pb_app_o)
      co_t <- .r4vn_score_calibration(yt, po_test)
      cs_a <- .r4vn_score_calibration(yb, pb_app_s)
      cs_t <- .r4vn_score_calibration(yt, ps_test)
      cutvals <- rep(NA_real_, 11L)
      names(cutvals) <- c("selected_cutoff",
                          "sens_app_cut", "sens_test_cut", "spec_app_cut", "spec_test_cut",
                          "ppv_app_cut", "ppv_test_cut", "npv_app_cut", "npv_test_cut",
                          "acc_app_cut", "acc_test_cut")
      if (!is.null(fitb$selected_cutoff) && is.finite(fitb$selected_cutoff)) {
        ca <- .r4vn_score_cut_metrics(yb, fitb$scores$clinical >= fitb$selected_cutoff)
        ct <- .r4vn_score_cut_metrics(yt, score_test >= fitb$selected_cutoff)
        getm_binary_boot <- function(tab, nm) tab$Estimate[match(nm, tab$Measure)]
        cutvals <- c(
          selected_cutoff = fitb$selected_cutoff,
          sens_app_cut = getm_binary_boot(ca, "Sensitivity"), sens_test_cut = getm_binary_boot(ct, "Sensitivity"),
          spec_app_cut = getm_binary_boot(ca, "Specificity"), spec_test_cut = getm_binary_boot(ct, "Specificity"),
          ppv_app_cut = getm_binary_boot(ca, "PPV"), ppv_test_cut = getm_binary_boot(ct, "PPV"),
          npv_app_cut = getm_binary_boot(ca, "NPV"), npv_test_cut = getm_binary_boot(ct, "NPV"),
          acc_app_cut = getm_binary_boot(ca, "Accuracy"), acc_test_cut = getm_binary_boot(ct, "Accuracy")
        )
      }
      out[[b]] <- c(
        auc_app_original = .r4vn_score_auc(yb, pb_app_o), auc_test_original = .r4vn_score_auc(yt, po_test),
        auc_app_score = .r4vn_score_auc(yb, pb_app_s), auc_test_score = .r4vn_score_auc(yt, ps_test),
        brier_app_original = mean((yb - pb_app_o)^2), brier_test_original = mean((yt - po_test)^2),
        brier_app_score = mean((yb - pb_app_s)^2), brier_test_score = mean((yt - ps_test)^2),
        calint_app_original = co_a[["intercept"]], calint_test_original = co_t[["intercept"]],
        calslope_app_original = co_a[["slope"]], calslope_test_original = co_t[["slope"]],
        calint_app_score = cs_a[["intercept"]], calint_test_score = cs_t[["intercept"]],
        calslope_app_score = cs_a[["slope"]], calslope_test_score = cs_t[["slope"]],
        cutvals
      )
    } else if (object$family == "cox") {
      mb_app_o <- as.numeric(stats::predict(fitb$models$original, type = "lp"))
      mb_app_s <- fitb$scores$clinical
      mt_o <- try(stats::predict(fitb, newdata = dat, type = "model_lp"), silent = TRUE)
      mt_s <- try(stats::predict(fitb, newdata = dat, type = "score"), silent = TRUE)
      if (inherits(mt_o, "try-error") || inherits(mt_s, "try-error")) next

      # Validate the time-dependent cutoff chosen inside each bootstrap sample.
      # The same bootstrap-derived numeric cutoff is then evaluated on the
      # original development sample. This preserves the full development
      # pipeline and therefore includes optimism from cutoff selection itself.
      cutvals <- rep(NA_real_, 11L)
      names(cutvals) <- c("selected_cutoff",
                          "sens_app_cut", "sens_test_cut", "spec_app_cut", "spec_test_cut",
                          "ppv_app_cut", "ppv_test_cut", "npv_app_cut", "npv_test_cut",
                          "acc_app_cut", "acc_test_cut")
      hz <- fitb$cutoff_time
      if (!is.null(fitb$selected_cutoff) && is.finite(fitb$selected_cutoff) &&
          !is.null(hz) && is.finite(hz)) {
        ca <- try(.r4vn_score_cut_metrics_cox(
          fitb$.development_data$.r4vn_score_time,
          fitb$.development_data$.r4vn_score_y,
          fitb$scores$clinical, fitb$selected_cutoff, hz
        ), silent = TRUE)
        ct <- try(.r4vn_score_cut_metrics_cox(
          object$.development_data$.r4vn_score_time,
          object$.development_data$.r4vn_score_y,
          mt_s, fitb$selected_cutoff, hz
        ), silent = TRUE)
        if (!inherits(ca, "try-error") && !inherits(ct, "try-error")) {
          getm_cox_boot <- function(tab, prefix) {
            jj <- grep(paste0("^", prefix), tab$Measure, ignore.case = TRUE)
            if (!length(jj)) return(NA_real_)
            as.numeric(tab$Estimate[jj[[1L]]])
          }
          cutvals <- c(
            selected_cutoff = fitb$selected_cutoff,
            sens_app_cut = getm_cox_boot(ca, "Sensitivity"), sens_test_cut = getm_cox_boot(ct, "Sensitivity"),
            spec_app_cut = getm_cox_boot(ca, "Specificity"), spec_test_cut = getm_cox_boot(ct, "Specificity"),
            ppv_app_cut = getm_cox_boot(ca, "PPV"), ppv_test_cut = getm_cox_boot(ct, "PPV"),
            npv_app_cut = getm_cox_boot(ca, "NPV"), npv_test_cut = getm_cox_boot(ct, "NPV"),
            acc_app_cut = getm_cox_boot(ca, "Accuracy"), acc_test_cut = getm_cox_boot(ct, "Accuracy")
          )
        }
      }

      out[[b]] <- c(
        c_app_original = .r4vn_score_cindex(fitb$.development_data$.r4vn_score_time, fitb$.development_data$.r4vn_score_y, mb_app_o)[1],
        c_test_original = .r4vn_score_cindex(object$.development_data$.r4vn_score_time, object$.development_data$.r4vn_score_y, mt_o)[1],
        c_app_score = .r4vn_score_cindex(fitb$.development_data$.r4vn_score_time, fitb$.development_data$.r4vn_score_y, mb_app_s)[1],
        c_test_score = .r4vn_score_cindex(object$.development_data$.r4vn_score_time, object$.development_data$.r4vn_score_y, mt_s)[1],
        cutvals
      )
    }
  }
  out <- out[!vapply(out, is.null, logical(1L))]
  if (!length(out)) return(data.frame(Note = "Bootstrap validation failed in all resamples."))
  M <- do.call(rbind, out)
  if (object$family == "logistic" || (object$family == "poisson" && object$binary_outcome)) {
    yy0 <- object$.development_data$.r4vn_score_y
    app_o_auc <- .r4vn_score_auc(yy0, object$predictions$original)
    app_s_auc <- .r4vn_score_auc(yy0, object$predictions$clinical)
    app_o_b <- mean((yy0 - object$predictions$original)^2)
    app_s_b <- mean((yy0 - object$predictions$clinical)^2)
    app_o_cal <- .r4vn_score_calibration(yy0, object$predictions$original)
    app_s_cal <- .r4vn_score_calibration(yy0, object$predictions$clinical)
    optimism <- c(
      AUC_original = mean(M[, "auc_app_original"] - M[, "auc_test_original"], na.rm = TRUE),
      Brier_original = mean(M[, "brier_app_original"] - M[, "brier_test_original"], na.rm = TRUE),
      CalInt_original = mean(M[, "calint_app_original"] - M[, "calint_test_original"], na.rm = TRUE),
      CalSlope_original = mean(M[, "calslope_app_original"] - M[, "calslope_test_original"], na.rm = TRUE),
      AUC_score = mean(M[, "auc_app_score"] - M[, "auc_test_score"], na.rm = TRUE),
      Brier_score = mean(M[, "brier_app_score"] - M[, "brier_test_score"], na.rm = TRUE),
      CalInt_score = mean(M[, "calint_app_score"] - M[, "calint_test_score"], na.rm = TRUE),
      CalSlope_score = mean(M[, "calslope_app_score"] - M[, "calslope_test_score"], na.rm = TRUE)
    )
    apparent <- c(app_o_auc, app_o_b, app_o_cal[["intercept"]], app_o_cal[["slope"]],
                  app_s_auc, app_s_b, app_s_cal[["intercept"]], app_s_cal[["slope"]])
    opt <- unname(optimism)
    vv <- data.frame(
      Model = rep(c("Original model", "Clinical score"), each = 4),
      Metric = rep(c("AUC", "Brier", "Calibration intercept", "Calibration slope"), 2),
      Apparent = apparent,
      Optimism = opt,
      Corrected = apparent - opt,
      Successful_bootstraps = nrow(M), stringsAsFactors = FALSE
    )
    if (!is.null(object$selected_cutoff) && is.finite(object$selected_cutoff) &&
        all(c("selected_cutoff", "sens_app_cut", "sens_test_cut") %in% colnames(M))) {
      cp <- object$tables$cutoff_performance_raw
      get_final_binary <- function(nm) cp$Estimate[match(nm, cp$Measure)]
      defs <- list(
        Sensitivity = c("sens_app_cut", "sens_test_cut"),
        Specificity = c("spec_app_cut", "spec_test_cut"),
        PPV = c("ppv_app_cut", "ppv_test_cut"),
        NPV = c("npv_app_cut", "npv_test_cut"),
        Accuracy = c("acc_app_cut", "acc_test_cut")
      )
      cr <- do.call(rbind, lapply(names(defs), function(nm) {
        cc <- defs[[nm]]
        op <- mean(M[, cc[[1L]]] - M[, cc[[2L]]], na.rm = TRUE)
        ap <- get_final_binary(nm)
        data.frame(Model = "Clinical score cutoff", Metric = nm, Apparent = ap,
                   Optimism = op, Corrected = ap - op, Successful_bootstraps = nrow(M),
                   stringsAsFactors = FALSE)
      }))
      vv <- rbind(vv, cr)
      bc <- M[, "selected_cutoff"]
      bc <- bc[is.finite(bc)]
      vv$Bootstrap_cutoff_median <- NA_real_
      vv$Bootstrap_cutoff_low <- NA_real_
      vv$Bootstrap_cutoff_high <- NA_real_
      if (length(bc)) {
        jj <- vv$Model == "Clinical score cutoff"
        q <- stats::quantile(bc, c(.025, .975), na.rm = TRUE, names = FALSE)
        vv$Bootstrap_cutoff_median[jj] <- stats::median(bc, na.rm = TRUE)
        vv$Bootstrap_cutoff_low[jj] <- q[[1L]]
        vv$Bootstrap_cutoff_high[jj] <- q[[2L]]
      }
    }
    vv
  } else {
    app_o <- .r4vn_score_cindex(object$.development_data$.r4vn_score_time, object$.development_data$.r4vn_score_y,
                                as.numeric(stats::predict(object$models$original, type = "lp")))[1]
    app_s <- .r4vn_score_cindex(object$.development_data$.r4vn_score_time, object$.development_data$.r4vn_score_y,
                                object$scores$clinical)[1]
    oo <- mean(M[, "c_app_original"] - M[, "c_test_original"], na.rm = TRUE)
    os <- mean(M[, "c_app_score"] - M[, "c_test_score"], na.rm = TRUE)
    vv <- data.frame(Model = c("Original model", "Clinical score"), Metric = "C-index",
                     Apparent = c(app_o, app_s), Optimism = c(oo, os), Corrected = c(app_o - oo, app_s - os),
                     Successful_bootstraps = nrow(M), stringsAsFactors = FALSE)

    if (!is.null(object$selected_cutoff) && is.finite(object$selected_cutoff) &&
        all(c("selected_cutoff", "sens_app_cut", "sens_test_cut") %in% colnames(M)) &&
        !is.null(object$tables$cutoff_performance_raw)) {
      cp <- object$tables$cutoff_performance_raw
      get_final_cox <- function(prefix) {
        jj <- grep(paste0("^", prefix), cp$Measure, ignore.case = TRUE)
        if (!length(jj)) return(NA_real_)
        as.numeric(cp$Estimate[jj[[1L]]])
      }
      defs <- list(
        `Sensitivity (IPCW)` = c("sens_app_cut", "sens_test_cut", "Sensitivity"),
        `Specificity (IPCW)` = c("spec_app_cut", "spec_test_cut", "Specificity"),
        `PPV (IPCW)` = c("ppv_app_cut", "ppv_test_cut", "PPV"),
        `NPV (IPCW)` = c("npv_app_cut", "npv_test_cut", "NPV"),
        `Accuracy (IPCW)` = c("acc_app_cut", "acc_test_cut", "Accuracy")
      )
      cr <- do.call(rbind, lapply(names(defs), function(nm) {
        cc <- defs[[nm]]
        op <- mean(M[, cc[[1L]]] - M[, cc[[2L]]], na.rm = TRUE)
        ap <- get_final_cox(cc[[3L]])
        data.frame(Model = "Clinical score cutoff", Metric = nm,
                   Apparent = ap, Optimism = op, Corrected = ap - op,
                   Successful_bootstraps = nrow(M), stringsAsFactors = FALSE)
      }))
      vv <- rbind(vv, cr)
      bc <- M[, "selected_cutoff"]
      bc <- bc[is.finite(bc)]
      vv$Bootstrap_cutoff_median <- NA_real_
      vv$Bootstrap_cutoff_low <- NA_real_
      vv$Bootstrap_cutoff_high <- NA_real_
      if (length(bc)) {
        jj <- vv$Model == "Clinical score cutoff"
        q <- stats::quantile(bc, c(.025, .975), na.rm = TRUE, names = FALSE)
        vv$Bootstrap_cutoff_median[jj] <- stats::median(bc, na.rm = TRUE)
        vv$Bootstrap_cutoff_low[jj] <- q[[1L]]
        vv$Bootstrap_cutoff_high[jj] <- q[[2L]]
      }
    }
    vv
  }
}

#' Build, simplify, validate and present a clinical/statistical scorecard
#'
#' @description
#' `tabscore()` converts a multivariable prediction model into a complete,
#' publication-ready scorecard. The function is designed for the full workflow,
#' not merely for rounding regression coefficients. In one call it can:
#'
#' * develop or accept a final prediction model;
#' * optionally select predictors from a candidate set;
#' * create score-ready categories for continuous predictors;
#' * derive an exact/model score and a simpler clinical integer score;
#' * produce a complete Predictor--Category--Point table and theoretical total
#'   score range;
#' * map every possible total score to predicted risk (or expected count/rate);
#' * select clinically/statistically useful score cutoffs;
#' * report diagnostic/prognostic properties at the selected cutoff;
#' * compare the original model, the model score and the clinical score;
#' * perform bootstrap internal validation of the *entire development pipeline*;
#' * optionally evaluate an external validation data set; and
#' * retain plot-ready data and prediction methods for deployment in R4VN Studio.
#'
#' Version 1 supports logistic regression, Cox proportional hazards regression,
#' and Poisson regression. Logistic and Cox models receive the most complete
#' discrimination/cutoff workflow. A Poisson model may be used for count/rate
#' scores; when its outcome is binary, ROC/cutoff summaries are also available.
#'
#' @param outcome Outcome variable. It can be an unquoted variable name, a
#'   character variable name, an already fitted `glm`/`coxph` model, or an R4VN
#'   regression result containing `raw$model` (for example from `logistic()` or
#'   `poisson()`). For Cox
#'   analysis this is the event/status indicator; supply follow-up time in
#'   `time=`. Binary outcomes may be numeric/logical/factor/character.
#' @param predictors Candidate predictors. Accepts a character vector, unquoted
#'   variables inside `c()`, or an R4VN `vars()` expression. Predictors are treated
#'   as model *variables* rather than individual dummy coefficients. When `vars()`
#'   is used, its R4VN declarations are honored: no prefix and `b2.`, `b3.`, ...
#'   declare categorical predictors and select the corresponding factor reference;
#'   `c.`, `q.` and `f.` declare numeric predictors as continuous for model fitting.
#'   Numeric variables carrying complete named value labels are also treated as
#'   categorical scorecard variables. Omit this argument when `outcome` is an
#'   already fitted model.
#' @param data Data frame. When omitted, `tabscore()` attempts to use the active
#'   R4VN data set created by `usedf()`. When `outcome` is an already fitted model,
#'   the stored model frame is the authoritative development sample; `data` is not
#'   used to refit or silently change that model.
#' @param family Model family: `auto`, `logistic`, `cox`, or `poisson`.
#'   With `auto`, the presence of `time=` selects Cox; otherwise a binary outcome
#'   selects logistic and a non-negative integer count selects Poisson.
#' @param event Event level for a binary outcome/status. For 0/1 outcomes the
#'   default is 1. For a two-level factor the default is its second factor level;
#'   for character outcomes it is the second observed non-missing value. Set this explicitly whenever the event
#'   direction matters, e.g. `event="Co"`.
#' @param time Cox follow-up-time variable, supplied as an unquoted or character
#'   variable name. This is the observed time, not the prediction horizon.
#' @param times Prediction horizons for Cox score-to-risk tables, e.g.
#'   `times=c(1,3,5)`. When omitted, useful event-time quantiles are selected.
#' @param cutoff_time Time horizon at which a Cox cutoff is evaluated. Defaults
#'   to the largest value in `times`, which is often the main clinical horizon.
#' @param select Predictor-selection strategy. `none` and `full` keep all
#'   supplied predictors. `backward`/`forward` use AIC stepwise selection.
#'   `purposeful` uses univariable screening, multivariable removal, confounding
#'   assessment and re-entry. `lasso` uses cross-validated `glmnet` at lambda.1se
#'   (falling back to lambda.min if necessary), then refits a standard model.
#' @param force Predictors that must remain in model-selection procedures. Accepts
#'   character names, `c(...)`, or `vars(...)`, using the same naming conventions
#'   as `predictors`.
#' @param exclude Candidate predictors to remove before model development. Accepts
#'   character names, `c(...)`, or `vars(...)`.
#' @param entry Univariable screening p-value for purposeful selection; default 0.25.
#' @param stay Multivariable retention/re-entry p-value for purposeful selection;
#'   default 0.10.
#' @param confound Relative coefficient-change threshold used to retain a variable
#'   as a confounder during purposeful selection; default 0.15 (15 percent).
#' @param cuts Optional named list of user-defined cut points for continuous
#'   predictors, e.g. `list(age=c(40,50,60), bmi=c(23,25,30))`. These are used in
#'   the score-ready model and therefore in the scorecard. As a convenience, a
#'   single character value `"easy"`, `"auto"`, `"quantile"`, or `"keep"`
#'   is accepted as an alias for `continuous=`. The original final model remains
#'   available for comparison.
#' @param continuous Handling of continuous predictors when the clinical scorecard
#'   is built. `easy` (default) uses quantile-informed cut points snapped to
#'   easy-to-use numbers; `auto` is an alias; `quantile` uses unsnapped empirical
#'   quantiles; `keep` leaves continuous predictors continuous in the score-ready
#'   model. A complete bedside Predictor--Category--Point table nevertheless needs
#'   explicit categories, so use `cuts=` when `keep` is requested.
#' @param bins Desired number of categories when automatic continuous-variable
#'   categorization is used. Default 4.
#' @param points Point-construction method. `auto` (or `clinical`/`integer`)
#'   searches for a small integer score that preserves score-ready model
#'   discrimination within `tolerance`. `model` or `pdo` uses model/PDO scaling.
#'   A named list or a data frame with columns Predictor, Category and Point may
#'   be supplied for a completely user-defined integer score. With the default
#'   `riskonly=TRUE`, manual points are shifted within each predictor so its
#'   minimum becomes zero; therefore even a supplied negative protective point is
#'   converted to an equivalent add-only risk score. With `riskonly=FALSE`,
#'   negative manual points are retained. `tabscore()` still calibrates, tests,
#'   compares and validates the manual score. When manual points refer to
#'   categories of a continuous predictor, also supply explicit `cuts=` so the
#'   category definitions remain fixed during validation. Named point vectors are
#'   matched to displayed category labels after normalizing common typographic
#'   equivalents, so ASCII input such as `40-49` and `>=60` also matches
#'   publication labels such as `40-49` and `>=60`. Unnamed vectors are matched in
#'   displayed category order.
#' @param pdo Points to double the effect on the model's log scale. For logistic
#'   regression this is Points to Double the Odds: `factor = pdo/log(2)`. For Cox
#'   the same number of points doubles hazard; for Poisson it doubles modeled rate.
#'   Default 20. The PDO/model score is retained separately from the clinical score.
#' @param maxscore Maximum preferred theoretical clinical score. `auto` searches
#'   compact totals (approximately 5--30 points). A numeric value constrains the
#'   theoretical range. This range is based on all possible scorecard categories,
#'   not merely the observed sample minimum/maximum.
#' @param simplify Logical. With `points="auto"`, TRUE (default) searches for
#'   a compact clinical integer score; FALSE uses rounded PDO/model points instead,
#'   which is useful when preserving model-scale resolution is more important than
#'   bedside compactness.
#' @param tolerance Maximum tolerated decrease in the primary discrimination
#'   metric when simplifying to the clinical score. For binary models this is AUC;
#'   for Cox it is C-index. The same threshold is also used to flag excessive loss
#'   caused by automatic predictor categorization. Default 0.01. If no compact score meets the tolerance,
#'   the best-performing candidate is retained and a warning is stored.
#' @param riskonly Logical. Default TRUE. Build the bedside score as a pure
#'   add-only risk score: each predictor is re-referenced to its lowest modeled
#'   risk category, all scoring effect ratios are at least 1, and all clinical
#'   points are non-negative. For example, if Male is the regression reference
#'   and Female has OR=0.50, the score representation becomes Female=0 points
#'   (score reference) and Male has risk-oriented OR=2.00 with positive points.
#'   This reparameterization does not change the fitted model, subject ranking,
#'   or predicted risks. Set FALSE only when a signed score with negative
#'   protective points is specifically desired.
#' @param scoreref Scoring reference. Default `"lowest"` chooses the lowest-risk
#'   category independently within every categorical predictor. `"model"` uses
#'   each fitted model reference and is mainly useful with `riskonly=FALSE`.
#'   A named list/vector can set explicit scoring references, for example
#'   `list(sex="Female", exercise="Yes")`. With `riskonly=TRUE`, an explicitly
#'   requested reference must be one of the predictor's lowest-risk categories;
#'   otherwise `tabscore()` stops because satisfying that reference would require
#'   negative risk points.
#' @param cutoff Selected cutoff principle: `iu`, `youden`, `risk`, `prevalence`,
#'   `sens`, `spec`, `cost`, `manual`, `refprob`, or `none`. All available methods
#'   are shown; this argument chooses the threshold used for final classification.
#' @param riskcut One or more clinically meaningful probability thresholds. With
#'   `cutoff="risk"`, the first value is converted to the nearest integer score.
#'   Multiple values create Low/Intermediate/High/etc. groups in the risk table.
#' @param cutoff_value User-specified integer score threshold for `cutoff="manual"`.
#' @param sens Minimum desired sensitivity for `cutoff="sens"`, e.g. 0.90.
#' @param spec Minimum desired specificity for `cutoff="spec"`, e.g. 0.90.
#' @param cost_fp Relative cost assigned to a false positive. Default 1.
#' @param cost_fn Relative cost assigned to a false negative. Example: `cost_fn=5`
#'   makes a false negative five times as costly as a false positive when
#'   `cost_fp=1`.
#' @param refprob Optional reference predicted probability for binary-outcome
#'   workflows: a numeric vector or a probability-variable name in `data`.
#'   `tabscore()` reports correlation, MAE, RMSE and mean difference between
#'   reference and score-derived probabilities. It is not used as a surrogate
#'   outcome and is not called a gold standard unless it truly is one.
#' @param refcut Probability threshold applied to `refprob`. With
#'   `cutoff="refprob"`, the score threshold that best reproduces this reference
#'   classification is selected, then evaluated against the actual outcome.
#' @param validate Internal validation: `bootstrap` (default) or `none`. Bootstrap
#'   validation reruns the whole development process inside each resample,
#'   including selection, automatic cuts and point simplification.
#' @param bootstrap Number of bootstrap resamples. Default 500. Use 1000 or more
#'   for a final analysis when feasible; small values are for code testing only.
#' @param validation Optional external validation data frame. The frozen final
#'   scorecard is applied without re-estimating cuts, points, calibration or the
#'   chosen score cutoff. Binary validation reports AUC, Brier, calibration and
#'   fixed-cutoff performance; Cox validation reports C-index, time-specific IPCW
#'   Brier and fixed-cutoff IPCW performance; count Poisson reports prediction error.
#' @param risktable Logical; create score-to-risk/score-to-expected-value table.
#' @param compare Logical; compare original model, model score and clinical score.
#' @param calibration Logical; prepare calibration data for plots.
#' @param decision Logical; prepare decision-curve net-benefit data for binary outcomes.
#' @param plot Logical. Default TRUE. Prepare all graphics supported by the chosen
#'   model family, embed every available graph directly in the HTML Viewer, and
#'   draw them into the interactive R/RStudio Plot history. Set FALSE when only
#'   tables are wanted. Plotting uses base R and does not add a required package.
#' @param show Logical. If TRUE, write a self-contained publication-oriented HTML
#'   result and open it in the RStudio Viewer (or the default browser).
#' @param console Logical. If TRUE, print a concise console summary.
#' @param seed Optional random seed used for LASSO, bootstrap fallback and validation. The default `NULL` does not set a seed.
#' @param ai Logical or endpoint name. If R4VN `aiask()` is available, request an
#'   optional AI interpretation after the statistical object is complete.
#'
#' @details
#' ## Missing data and development sample
#'
#' When `tabscore()` develops a model from candidate predictors, it uses one
#' complete-case development sample across the outcome/status, Cox time (when
#' applicable), and all candidate predictors supplied before model selection.
#' This keeps candidate models comparable but can reduce sample size when many
#' predictors have missing values. Perform the intended imputation or missing-data
#' strategy before `tabscore()` when complete-case analysis is inappropriate.
#'
#' ## R4VN vars() declarations and reference categories
#'
#' `tabscore()` understands the R4VN variable culture rather than merely stripping
#' prefixes. For example, `vars(c.age, b2.sex, smoking)` fits age continuously,
#' treats sex as categorical with its second factor level as model reference, and
#' treats smoking as categorical with its first factor level as reference. This
#' affects the fitted model table and model-selection calculations. Point assignment
#' itself is then shifted within each predictor so the lowest-risk category receives
#' zero automatic points; consequently the zero-point category does not have to be
#' the regression reference category. The final prediction-model table preserves
#' the original statistical reference and may therefore legitimately show OR/HR/RR
#' below 1. The separate `risk_orientation` table shows the scoring contrast after
#' re-referencing; with `riskonly=TRUE`, every displayed scoring ratio is >=1 and
#' the clinical score contains only zero or positive points. With ordinary
#' `c(age, sex)` syntax, data type is inferred from the columns instead of imposing
#' R4VN categorical declarations.
#'
#' ## Original model, model score and clinical score
#'
#' `tabscore()` distinguishes three objects. The **original model** is the selected
#' or fixed model using the original predictor representation. The **model score**
#' is a monotone point transformation of the score-ready model linear predictor.
#' With `pdo=20`, a 20-point increase doubles odds (logistic), hazard (Cox), or
#' modeled rate (Poisson). The **clinical score** uses small integer points and is
#' the score shown in the clean Predictor--Category--Point table.
#'
#' Intercepts are never artificially divided among predictors; they remain in the
#' risk mapping. Within each predictor, the lowest modeled contribution is shifted
#' to zero before automatic non-negative points are assigned. Thus a zero-point
#' category need not be the regression reference if another category has lower risk.
#'
#' ## Protective factors and add-only risk scoring
#'
#' With `riskonly=TRUE`, categorical contributions are transformed predictor by
#' predictor as `beta_score = beta_category - min(beta_categories)`. Therefore
#' `exp(beta_score) >= 1`. For a binary predictor with an original protective
#' contrast OR=0.50, reversing the scoring contrast gives 1/0.50=2.00 for the
#' higher-risk category. For multi-level predictors the same minimum-risk rebasing
#' is used; the procedure is not `abs(beta)`, which can distort category ordering.
#' Cox HR and Poisson RR/IRR are handled identically on their log-effect scales.
#' For a continuous coefficient retained without categories, a negative coefficient
#' is described technically as risk per unit decrease; a complete bedside integer
#' score still requires explicit/automatic categories. The transformation changes
#' only the score origin/reference: the original fitted model and its absolute-risk
#' predictions remain untouched.
#'
#' ## Starting from an already fitted model
#'
#' A base `glm`/`coxph` model or an R4VN regression result containing `raw$model`
#' may be passed as `outcome`. In that workflow the fitted predictor set is fixed,
#' so `select`, `force` and `exclude` are not used. This version deliberately
#' rejects interactions, transformed/spline terms, no-intercept models, non-unit
#' analysis weights, and non-zero offsets/exposures rather than silently changing
#' the fitted model during score simplification. Represent required transformed
#' predictors as explicit columns and refit before calling `tabscore()`, or provide
#' a manual point system.
#'
#' ## Automatic and manual cut points
#'
#' Automatic scorecard categorization is a simplification step rather than part of
#' the original continuous model. `continuous="easy"` uses empirical quantiles and
#' snaps thresholds toward simple numbers. Prefer clinically established `cuts=`
#' when available. Because automatic cuts are data-driven, bootstrap validation
#' repeats the cut-selection step inside each resample.
#'
#' ## Theoretical total score
#'
#' The displayed range is calculated from all category combinations implied by the
#' scorecard, not from the smallest/largest observed subject score. This keeps the
#' bedside score stable in new data.
#'
#' ## Score-to-risk conversion
#'
#' For logistic models, the clinical score is recalibrated by
#' `logit(P)=a+b*Score`; every possible total score receives a predicted probability
#' and 95 percent CI. This recalibration is important because categorization and
#' integer rounding mean the clinical score is no longer exactly the original LP.
#' Cox scorecards use a one-predictor Cox calibration model to provide risk at each
#' `times=` horizon. Poisson count scores provide expected count/rate and CI.
#'
#' ## Cutoff methods
#'
#' * `youden`: maximize sensitivity + specificity - 1.
#' * `iu`: minimize `abs(sensitivity-AUC)+abs(specificity-AUC)`.
#' * `risk`: map a clinical probability in `riskcut` to an integer score.
#' * `prevalence`: use development event prevalence as probability threshold,
#'   reproducing a common legacy workflow; it is not automatically clinically best.
#' * `sens`: among thresholds meeting `sens`, maximize specificity.
#' * `spec`: among thresholds meeting `spec`, maximize sensitivity.
#' * `cost`: minimize `cost_fn*FN + cost_fp*FP`.
#' * `manual`: use `cutoff_value`.
#' * `refprob`: use `refprob` and `refcut` to reproduce a reference probability rule.
#'
#' Cox support in this version is for standard right-censored proportional-hazards
#' models. Start-stop/time-dependent, multi-state, competing-risk and other complex
#' survival structures are not silently reduced to a simple integer score.
#' Cox ROC/cutoff calculations at `cutoff_time` use cumulative/dynamic IPCW
#' sensitivity and specificity with the censoring distribution estimated by
#' Kaplan-Meier. IPCW cutoff sensitivity, specificity, PPV, NPV, accuracy and
#' likelihood ratios are point estimates in the apparent cutoff table; unlike the
#' ordinary binary-outcome table, simple binomial confidence intervals are not
#' reported because censoring weights make them inappropriate. Full-pipeline
#' bootstrap validation supplies optimism-corrected cutoff performance and the
#' empirical stability interval of the selected score threshold. If an intervention
#' probability is known, `riskcut` is generally easier to interpret clinically than
#' a purely statistical cutoff.
#'
#' ## Model versus score comparison
#'
#' A score is not accepted merely because an AUC difference is non-significant.
#' Binary comparisons include AUC with 95 percent CI, Brier score, calibration
#' intercept/slope, and paired AUC difference. Cox comparisons include C-index and
#' time-specific IPCW Brier scores. Decision-curve data compare net benefit across
#' probability thresholds.
#'
#' ## Bootstrap internal validation
#'
#' Each bootstrap resample repeats predictor selection, automatic categorization,
#' point derivation, and the requested cutoff-selection rule. Apparent performance
#' is compared with performance when the bootstrap-derived score and cutoff are
#' applied to the original sample. For binary outcomes the validation table includes
#' optimism-corrected AUC, Brier score, calibration intercept/slope, sensitivity,
#' specificity, PPV, NPV, accuracy, and bootstrap cutoff stability. For Cox models
#' the validation table includes optimism-corrected C-index and, when a cutoff is
#' requested, time-dependent IPCW sensitivity, specificity, PPV, NPV and accuracy
#' at `cutoff_time` plus bootstrap cutoff stability. Mean optimism is subtracted
#' from the apparent final performance. This is more rigorous than bootstrapping
#' a fixed already-developed score.
#'
#' ## External validation and deployment
#'
#' Supply `validation=` to apply the *frozen* developed score to a separate data
#' set. Predictor cut points, integer points, risk mapping, and the selected score
#' threshold are not re-optimized in the external data. This avoids turning an
#' external validation into a second development exercise. After development,
#' `predict(score_object, newdata, type="all")` returns bedside score, predicted
#' risk/value and risk group where applicable. New factor values that were absent
#' from the development scorecard cannot be scored and therefore yield missing
#' score/risk rather than being silently assigned zero points.
#'
#' ## Viewer, plots and package requirements
#'
#' With `show=TRUE`, the Viewer contains the publication tables and, when
#' `plot=TRUE`, every graph that is actually available for the fitted family.
#' Logistic/binary scorecards can show score-to-risk, ROC, calibration, decision
#' curve and score-distribution plots. Cox scorecards show score-to-risk across
#' requested horizons, time-dependent ROC at `cutoff_time`, and score
#' distribution. Poisson count scorecards show score-to-expected-value and the
#' observed-count distribution. The Viewer plots are rendered with base R, using
#' a system sans-serif font and an embedded raster image when possible; this keeps
#' the report self-contained and avoids a ggplot2/htmlwidgets dependency. Use
#' `plot(result)` to draw all available graphs in the RStudio Plots pane, or
#' `plot(result, which="roc")`, for example, to draw one.
#'
#' The core logistic and Poisson workflows use base/recommended R only.
#' `survival` is needed only for Cox scorecards, `glmnet` only for
#' `select="lasso"`, and `pROC` is optional because R4VN has a base-R fallback
#' for AUC calculations and paired AUC comparison. Thus users do not need to
#' install a large collection of packages for ordinary `tabscore()` analyses.
#'
#' ## Modeling cautions
#'
#' Prediction modeling is not equivalent to retaining only p<0.05 predictors. Use
#' subject-matter knowledge and `force=` for essential variables. Data-driven cuts
#' and cutoffs can overfit and should be validated. Interactions, spline bases,
#' time-varying Cox effects, competing risks and machine-learning distillation are
#' not silently converted into a bedside integer score in this version.
#'
#' @return An object of class `r4vn_tabscore` with fitted models, score rules,
#' development scores/predictions, selected cutoff, theoretical range, publication
#' tables, technical tables, validation results and plot-ready data. Important
#' elements include `models`, `scores`, `tables`, `publication_tables`,
#' `selected_predictors`, `selected_cutoff`, `score_range`, `riskonly`,
#' `effect_measure`, `plots`, `plot_titles` and `plot_data`.
#' `tables$risk_orientation` explicitly compares
#' the score-ready model-reference effect ratio with the risk-oriented scoring ratio. The
#' `publication_tables` list contains only ready-to-export non-NULL tables and can
#' be passed directly to `tabexport()`. `predict()` can then score new
#' patients without re-estimating the scorecard.
#'
#' @examples
#' set.seed(2026)
#' n <- 220
#' d <- data.frame(
#'   age = round(rnorm(n, 52, 12)),
#'   bmi = round(rnorm(n, 24, 4), 1),
#'   hypertension = factor(rbinom(n, 1, .30), 0:1, c("No", "Yes")),
#'   smoking = factor(rbinom(n, 1, .25), 0:1, c("No", "Yes")),
#'   alcohol = factor(rbinom(n, 1, .20), 0:1, c("No", "Yes"))
#' )
#' lp <- -4.2 + .04*d$age + .06*(d$bmi - 24) +
#'   .8*(d$hypertension == "Yes") + .6*(d$smoking == "Yes")
#' d$event <- rbinom(n, 1, plogis(lp))
#'
#' # 1. Simplest publication-ready logistic scorecard.
#' # show=TRUE and plot=TRUE are the user-facing defaults.
#' s1 <- tabscore(
#'   event, c(age, bmi, hypertension, smoking), data=d,
#'   validate="none", show=FALSE, plot=FALSE
#' )
#' s1$tables$scorecard
#' s1$tables$risk
#' s1$tables$comparison
#'
#' # 2. R4VN variable declarations: continuous variables and chosen references.
#' s2 <- tabscore(
#'   event, vars(c.age, c.bmi, b2.hypertension, b2.smoking), data=d,
#'   validate="none", show=FALSE, plot=FALSE
#' )
#' s2$tables$model
#' s2$tables$risk_orientation
#'
#' # 3. Clinically prespecified cut points.
#' s3 <- tabscore(
#'   event, c(age, bmi, hypertension, smoking), data=d,
#'   cuts=list(age=c(40,50,60), bmi=c(23,25,30)),
#'   validate="none", show=FALSE, plot=FALSE
#' )
#'
#' # 4. Apply the frozen scorecard to new patients.
#' newp <- data.frame(
#'   age=c(45,68), bmi=c(24,29),
#'   hypertension=factor(c("No","Yes"), levels=c("No","Yes")),
#'   smoking=factor(c("Yes","No"), levels=c("No","Yes"))
#' )
#' predict(s3, newp, type="all")
#'
#' \donttest{
#' # 5. Purposeful selection; force variables that must remain clinically.
#' s5 <- tabscore(
#'   event, c(age,bmi,hypertension,smoking,alcohol), data=d,
#'   select="purposeful", force=c("age","hypertension"),
#'   validate="none", show=FALSE, plot=FALSE
#' )
#'
#' # 6. Backward or forward AIC selection.
#' s6a <- tabscore(event, c(age,bmi,hypertension,smoking,alcohol), data=d,
#'                 select="backward", validate="none", show=FALSE, plot=FALSE)
#' s6b <- tabscore(event, c(age,bmi,hypertension,smoking,alcohol), data=d,
#'                 select="forward", validate="none", show=FALSE, plot=FALSE)
#'
#' # 7. LASSO is optional and only needs glmnet for this selection method.
#' if (requireNamespace("glmnet", quietly=TRUE)) {
#'   s7 <- tabscore(event, c(age,bmi,hypertension,smoking,alcohol), data=d,
#'                  select="lasso", validate="none", show=FALSE, plot=FALSE)
#' }
#'
#' # 8. Compact score versus PDO/model-scale points.
#' s8a <- tabscore(event, c(age,hypertension,smoking), data=d,
#'                 cuts=list(age=c(40,50,60)), maxscore=10,
#'                 validate="none", show=FALSE, plot=FALSE)
#' s8b <- tabscore(event, c(age,hypertension,smoking), data=d,
#'                 cuts=list(age=c(40,50,60)), points="pdo", pdo=20,
#'                 validate="none", show=FALSE, plot=FALSE)
#'
#' # 9. Completely manual bedside points; R4VN still calibrates and validates it.
#' s9 <- tabscore(
#'   event, c(age,hypertension,smoking), data=d,
#'   cuts=list(age=c(40,50,60)),
#'   points=list(
#'     age=c("<40"=0, "40-49"=1, "50-59"=2, ">=60"=3),
#'     hypertension=c("No"=0,"Yes"=2),
#'     smoking=c("No"=0,"Yes"=1)
#'   ), validate="none", show=FALSE, plot=FALSE
#' )
#'
#' # 10. Common cutoff rules.
#' s10_iu <- tabscore(event, c(age,hypertension,smoking), data=d,
#'                    cutoff="iu", validate="none", show=FALSE, plot=FALSE)
#' s10_youden <- tabscore(event, c(age,hypertension,smoking), data=d,
#'                        cutoff="youden", validate="none", show=FALSE, plot=FALSE)
#' s10_sens <- tabscore(event, c(age,hypertension,smoking), data=d,
#'                      cutoff="sens", sens=.90, validate="none", show=FALSE, plot=FALSE)
#' s10_cost <- tabscore(event, c(age,hypertension,smoking), data=d,
#'                      cutoff="cost", cost_fn=5, cost_fp=1,
#'                      validate="none", show=FALSE, plot=FALSE)
#'
#' # 11. Clinically meaningful probability threshold and multiple risk groups.
#' s11 <- tabscore(event, c(age,hypertension,smoking), data=d,
#'                 riskcut=c(.05,.10,.20), cutoff="risk",
#'                 validate="none", show=FALSE, plot=FALSE)
#' s11$tables$risk
#'
#' # 12. Compare the score with an existing/reference probability.
#' d$reference_risk <- plogis(-4 + .04*d$age + .7*(d$hypertension == "Yes"))
#' s12 <- tabscore(event, c(age,hypertension,smoking), data=d,
#'                 refprob=reference_risk, refcut=.10, cutoff="refprob",
#'                 validate="none", show=FALSE, plot=FALSE)
#' s12$tables$reference_probability
#'
#' # 13. Convert an already fitted logistic model.
#' m13 <- glm(event ~ age + hypertension + smoking, data=d, family=binomial())
#' s13 <- tabscore(m13, validate="none", show=FALSE, plot=FALSE)
#'
#' # 14. External validation with a frozen scorecard.
#' dev <- d[1:150, ]
#' val <- d[151:nrow(d), ]
#' s14 <- tabscore(event, c(age,hypertension,smoking), data=dev,
#'                 cuts=list(age=c(40,50,60)), validation=val,
#'                 validate="none", show=FALSE, plot=FALSE)
#' s14$tables$external_validation
#'
#' # 15. Full-pipeline bootstrap validation. B=20 is only a quick code check;
#' # use bootstrap=500 or more for the final report.
#' s15 <- tabscore(event, c(age,hypertension,smoking), data=d,
#'                 cuts=list(age=c(40,50,60)),
#'                 validate="bootstrap", bootstrap=20,
#'                 show=FALSE, plot=FALSE)
#' s15$tables$validation
#'
#' # 16. Cox scorecard: survival is the only package required for this family.
#' if (requireNamespace("survival", quietly=TRUE)) {
#'   ds <- d
#'   true_t <- rexp(nrow(ds), rate=exp(-3 + .02*ds$age +
#'                   .6*(ds$hypertension == "Yes")))
#'   censor_t <- rexp(nrow(ds), rate=.08)
#'   ds$status <- as.integer(true_t <= censor_t)
#'   ds$ftime <- pmin(true_t, censor_t)
#'   sc <- tabscore(status, c(age,hypertension,smoking), data=ds,
#'                  family="cox", time=ftime, times=c(1,3,5), cutoff_time=5,
#'                  cuts=list(age=c(40,50,60)),
#'                  validate="none", show=FALSE, plot=FALSE)
#'   sc$tables$risk
#'   sc$tables$time_brier
#' }
#'
#' # 17. Poisson count scorecard.
#' dp <- d
#' dp$count <- rpois(nrow(dp), exp(-1 + .015*dp$age + .35*(dp$smoking == "Yes")))
#' sp <- tabscore(count, c(age,smoking), data=dp, family="poisson",
#'                cuts=list(age=c(40,50,60)), validate="none",
#'                show=FALSE, plot=FALSE)
#' sp$tables$risk
#'
#' # 18. Every available plot; Viewer includes the same figures when show=TRUE.
#' sv <- tabscore(event, c(age,hypertension,smoking), data=d,
#'                cuts=list(age=c(40,50,60)), validate="none",
#'                show=FALSE, plot=TRUE)
#' plot(sv, which="risk")
#' plot(sv, which="roc")
#' plot(sv, which="calibration")
#' plot(sv, which="decision")
#' plot(sv, which="distribution")
#' plot(sv)  # all available plots in Plot history
#'
#' # 19. Protective predictors: default risk-only coding versus a signed score.
#' dr <- d
#' dr$exercise <- factor(rbinom(nrow(dr), 1, .55), 0:1, c("No", "Yes"))
#' dr$event2 <- rbinom(nrow(dr), 1,
#'                     plogis(-2.5 + .04*dr$age - .8*(dr$exercise == "Yes")))
#' srisk <- tabscore(event2, c(age,exercise), data=dr,
#'                   cuts=list(age=c(40,50,60)), riskonly=TRUE,
#'                   validate="none", show=FALSE, plot=FALSE)
#' ssigned <- tabscore(event2, c(exercise), data=dr,
#'                     points=list(exercise=c("No"=0,"Yes"=-2)),
#'                     riskonly=FALSE, scoreref="model",
#'                     validate="none", show=FALSE, plot=FALSE)
#'
#' # 20. Additional cutoff strategies.
#' s20_prev <- tabscore(event, c(age,hypertension,smoking), data=d,
#'                      cutoff="prevalence", validate="none",
#'                      show=FALSE, plot=FALSE)
#' s20_spec <- tabscore(event, c(age,hypertension,smoking), data=d,
#'                      cutoff="spec", spec=.90, validate="none",
#'                      show=FALSE, plot=FALSE)
#' s20_manual <- tabscore(event, c(age,hypertension,smoking), data=d,
#'                        cutoff="manual", cutoff_value=3, validate="none",
#'                        show=FALSE, plot=FALSE)
#' s20_none <- tabscore(event, c(age,hypertension,smoking), data=d,
#'                      cutoff="none", validate="none",
#'                      show=FALSE, plot=FALSE)
#'
#' # 21. Quantile-based automatic categorization.
#' s21 <- tabscore(event, c(age,bmi,hypertension,smoking), data=d,
#'                 continuous="quantile", bins=4,
#'                 validate="none", show=FALSE, plot=FALSE)
#'
#' # 22. An R4VN logistic() result can be converted directly as well.
#' rfit <- logistic(event, vars=vars(c.age, hypertension, smoking),
#'                  data=d, show=FALSE)
#' s22 <- tabscore(rfit, validate="none", show=FALSE, plot=FALSE)
#'
#' # 23. Prediction outputs after the scorecard is frozen.
#' predict(s3, newp, type="score")
#' predict(s3, newp, type="risk")
#' predict(s3, newp, type="group")
#' predict(s3, newp, type="model")
#'
#' # 24. Export all publication-ready tables without another tabscore-specific
#' # dependency. Word/Excel writers are only needed when those formats are chosen.
#' # tabexport(sv$publication_tables, export=c("html","docx","xlsx"),
#' #           file="tabscore_report", open=FALSE)
#' }
#' @export
#' @family prediction models
#' @family scorecards
#' @family R4VN tables
#' @seealso `predict.r4vn_tabscore`, `plot.r4vn_tabscore`
tabscore <- function(
    outcome,
    predictors = NULL,
    data = NULL,
    family = c("auto", "logistic", "cox", "poisson"),
    event = NULL,
    time = NULL,
    times = NULL,
    cutoff_time = NULL,
    select = c("none", "full", "backward", "forward", "purposeful", "lasso"),
    force = NULL,
    exclude = NULL,
    entry = 0.25,
    stay = 0.10,
    confound = 0.15,
    cuts = NULL,
    continuous = c("easy", "auto", "quantile", "keep"),
    bins = 4L,
    points = "auto",
    pdo = 20,
    maxscore = "auto",
    simplify = TRUE,
    tolerance = 0.01,
    riskonly = TRUE,
    scoreref = "lowest",
    cutoff = c("iu", "youden", "risk", "prevalence", "sens", "spec", "cost", "manual", "refprob", "none"),
    riskcut = NULL,
    cutoff_value = NULL,
    sens = NULL,
    spec = NULL,
    cost_fp = 1,
    cost_fn = 1,
    refprob = NULL,
    refcut = NULL,
    validate = c("bootstrap", "none"),
    bootstrap = 500L,
    validation = NULL,
    risktable = TRUE,
    compare = TRUE,
    calibration = TRUE,
    decision = TRUE,
    plot = TRUE,
    show = TRUE,
    console = FALSE,
    seed = NULL,
    ai = FALSE
) {
  env <- parent.frame()
  family <- match.arg(family)
  select <- match.arg(select)
  continuous <- match.arg(continuous)
  if (!is.logical(riskonly) || length(riskonly) != 1L || is.na(riskonly))
    stop("`riskonly` must be TRUE or FALSE.", call. = FALSE)
  .r4vn_score_validate_scoreref(scoreref)
  if (isTRUE(riskonly) && is.character(scoreref) && length(scoreref) == 1L &&
      is.null(names(scoreref)) && identical(tolower(scoreref), "model")) {
    warning("`riskonly=TRUE, scoreref='model'` can require negative points when the model reference is protective. `scoreref='lowest'` is used instead.", call. = FALSE)
    scoreref <- "lowest"
  }
  if (is.character(cuts) && length(cuts) == 1L) {
    cm <- match.arg(tolower(cuts), c("easy", "auto", "quantile", "keep"))
    continuous <- cm
    cuts <- NULL
  }
  cutoff <- match.arg(cutoff)
  validate <- match.arg(validate)
  if (!is.numeric(pdo) || length(pdo) != 1L || !is.finite(pdo) || pdo <= 0)
    stop("pdo must be a positive number.", call. = FALSE)
  if (!is.numeric(tolerance) || length(tolerance) != 1L || !is.finite(tolerance) || tolerance < 0)
    stop("tolerance must be one non-negative finite number.", call. = FALSE)
  if (!is.numeric(bins) || length(bins) != 1L || !is.finite(bins) || bins < 2)
    stop("bins must be an integer-like value of at least 2.", call. = FALSE)
  bins <- as.integer(round(bins))
  for (zz in c(entry, stay)) if (!is.numeric(zz) || length(zz) != 1L || !is.finite(zz) || zz < 0 || zz > 1)
    stop("entry and stay must be probabilities between 0 and 1.", call. = FALSE)
  if (!is.numeric(confound) || length(confound) != 1L || !is.finite(confound) || confound < 0)
    stop("confound must be a non-negative finite relative-change threshold.", call. = FALSE)
  if (!is.numeric(cost_fp) || length(cost_fp) != 1L || !is.finite(cost_fp) || cost_fp < 0 ||
      !is.numeric(cost_fn) || length(cost_fn) != 1L || !is.finite(cost_fn) || cost_fn < 0 ||
      (cost_fp == 0 && cost_fn == 0))
    stop("cost_fp and cost_fn must be non-negative finite numbers and cannot both be zero.", call. = FALSE)
  if (!is.null(riskcut) && (!is.numeric(riskcut) || any(!is.finite(riskcut)) || any(riskcut <= 0 | riskcut >= 1)))
    stop("riskcut must contain probability thresholds strictly between 0 and 1.", call. = FALSE)
  if (!is.null(sens) && (!is.numeric(sens) || length(sens) != 1L || !is.finite(sens) || sens <= 0 || sens > 1))
    stop("sens must be a probability in (0, 1].", call. = FALSE)
  if (!is.null(spec) && (!is.numeric(spec) || length(spec) != 1L || !is.finite(spec) || spec <= 0 || spec > 1))
    stop("spec must be a probability in (0, 1].", call. = FALSE)
  if (!is.null(refcut) && (!is.numeric(refcut) || length(refcut) != 1L || !is.finite(refcut) || refcut <= 0 || refcut >= 1))
    stop("refcut must be a probability strictly between 0 and 1.", call. = FALSE)
  if (!is.null(cutoff_value) && (!is.numeric(cutoff_value) || length(cutoff_value) != 1L || !is.finite(cutoff_value)))
    stop("cutoff_value must be one finite score value.", call. = FALSE)
  if (!is.null(times) && (!is.numeric(times) || any(!is.finite(times)) || any(times <= 0)))
    stop("times must contain positive finite Cox prediction horizons.", call. = FALSE)
  if (!is.null(cutoff_time) && (!is.numeric(cutoff_time) || length(cutoff_time) != 1L || !is.finite(cutoff_time) || cutoff_time <= 0))
    stop("cutoff_time must be one positive finite Cox horizon.", call. = FALSE)
  if (!identical(maxscore, "auto") && (!is.numeric(maxscore) || length(maxscore) != 1L || !is.finite(maxscore) || maxscore < 2))
    stop("maxscore must be 'auto' or one finite value of at least 2.", call. = FALSE)
  if (!is.numeric(bootstrap) || length(bootstrap) != 1L || !is.finite(bootstrap) || bootstrap < 0)
    stop("bootstrap must be a non-negative integer-like number.", call. = FALSE)
  bootstrap <- as.integer(round(bootstrap))
  if (cutoff == "risk" && (is.null(riskcut) || !length(riskcut)))
    stop("cutoff='risk' requires riskcut=, for example riskcut=0.10.", call. = FALSE)
  if (cutoff == "sens" && is.null(sens)) stop("cutoff='sens' requires sens=.", call. = FALSE)
  if (cutoff == "spec" && is.null(spec)) stop("cutoff='spec' requires spec=.", call. = FALSE)
  if (cutoff == "manual" && is.null(cutoff_value)) stop("cutoff='manual' requires cutoff_value=.", call. = FALSE)
  out_expr <- substitute(outcome)
  out_val <- try(eval(out_expr, env), silent = TRUE)
  supplied_fit <- if (inherits(out_val, "try-error")) NULL else .r4vn_score_extract_model(out_val)
  supplied_model <- !is.null(supplied_fit)
  source_event_label <- if (!inherits(out_val, "try-error") && is.list(out_val) &&
                             !is.null(out_val$raw) && is.list(out_val$raw)) out_val$raw$event else NULL
  source_vcov <- if (!inherits(out_val, "try-error") && is.list(out_val) &&
                       !is.null(out_val$raw) && is.list(out_val$raw)) out_val$raw$vcov else NULL
  refprob_expr <- substitute(refprob)
  time_expr <- substitute(time)
  force_expr <- substitute(force)
  exclude_expr <- substitute(exclude)
  time_supplied <- !identical(time_expr, quote(NULL))
  refprob_supplied <- !identical(refprob_expr, quote(NULL))
  if (cutoff == "refprob" && (!refprob_supplied || is.null(refcut)))
    stop("cutoff='refprob' requires both refprob= and refcut=.", call. = FALSE)
  force_names <- if (identical(force_expr, quote(NULL))) character() else .r4vn_score_predictor_names(force_expr, env)
  exclude_names <- if (identical(exclude_expr, quote(NULL))) character() else .r4vn_score_predictor_names(exclude_expr, env)

  source_model <- NULL
  outcome_name <- NULL
  time_name <- NULL
  event_info <- NULL
  original_row_index <- NULL
  refprob_full <- NULL
  refprob_name <- NULL
  data_raw_model <- NULL
  predictor_spec <- NULL

  if (supplied_model) {
    source_model <- supplied_fit
    if (!identical(select, "none") || length(force_names) || length(exclude_names)) {
      stop("When outcome is an already fitted model, its predictor set is fixed. Do not use select=, force=, or exclude=; refit the model first if the predictor set must change.", call. = FALSE)
    }
    tt0 <- stats::terms(source_model)
    if (isFALSE(attr(tt0, "intercept") == 1L))
      stop("Fitted models without an intercept are not converted automatically in this version.", call. = FALSE)
    term_labels0 <- attr(tt0, "term.labels")
    unsupported_terms <- term_labels0[grepl(":|\\(|\\)|\\^|\\*", term_labels0)]
    if (length(unsupported_terms)) {
      stop(
        "A fitted-model scorecard currently accepts simple main-effect terms only. ",
        "Interactions, transformations, polynomial/spline terms, offsets/strata, and other constructed terms ",
        "must first be represented as explicit data columns or the score must be specified manually. ",
        "Unsupported term(s): ", paste(unsupported_terms, collapse = ", "),
        call. = FALSE
      )
    }
    mf <- stats::model.frame(source_model)
    mw <- stats::model.weights(mf)
    if (!is.null(mw) && any(is.finite(mw) & abs(mw - 1) > 1e-10))
      stop("A fitted model with non-unit analysis weights cannot yet be automatically simplified without changing its estimation. Refit tabscore() from the original outcome/predictors after creating the intended score-ready variables, or use a manual score.", call. = FALSE)
    mo <- stats::model.offset(mf)
    if (!is.null(mo) && any(is.finite(mo) & abs(mo) > 1e-10))
      stop("A fitted model containing a non-zero offset/exposure cannot yet be automatically converted to the complete score-to-risk table in this version.", call. = FALSE)
    if (!nrow(mf)) stop("The supplied fitted model does not contain a usable model frame.", call. = FALSE)
    data_raw <- as.data.frame(mf)
    response <- stats::model.response(mf)
    lhs0 <- stats::formula(source_model)[[2L]]
    if (is.symbol(lhs0)) outcome_name <- as.character(lhs0)
    lhs_vars0 <- all.vars(lhs0)
    predictors0 <- unique(all.vars(stats::delete.response(stats::terms(source_model))))
    predictors0 <- predictors0[predictors0 %in% names(data_raw)]
    if (!length(predictors0)) predictors0 <- names(mf)[-1L]
    predictors0 <- setdiff(predictors0, exclude_names)
    labels <- setNames(vapply(predictors0, function(v)
      .r4vn_score_label(data_raw[[v]], v), character(1L)), predictors0)

    if (inherits(source_model, "coxph")) {
      family <- "cox"
      if (!requireNamespace("survival", quietly = TRUE))
        stop("Package 'survival' is required.", call. = FALSE)
      if (!inherits(response, "Surv"))
        stop("The coxph model response is not a Surv object.", call. = FALSE)
      stype <- attr(response, "type")
      if (!is.null(stype) && !stype %in% c("right"))
        stop("This version converts standard right-censored Cox models only; start-stop, interval, multi-state, or other Surv types require a specialized score workflow.", call. = FALSE)
      if (length(lhs_vars0) >= 2L) {
        time_name <- lhs_vars0[[1L]]
        outcome_name <- lhs_vars0[[length(lhs_vars0)]]
      }
      data_raw$.r4vn_score_time <- as.numeric(response[, 1L])
      data_raw$.r4vn_score_y <- as.integer(response[, ncol(response)])
      event_info <- list(event = 1, levels = c(0, 1))
    } else {
      famname <- source_model$family$family
      if (famname == "binomial") family <- "logistic"
      else if (famname == "poisson") family <- "poisson"
      else stop("Only binomial/logistic and Poisson glm models are supported.", call. = FALSE)
      if (family == "logistic") {
        event_fit <- if (is.numeric(response) || is.logical(response)) NULL else event
        bi <- .r4vn_score_binary(response, event_fit, "model outcome")
        data_raw$.r4vn_score_y <- bi$y
        event_info <- bi
        if (!is.null(source_event_label) && length(source_event_label)) event_info$event <- as.character(source_event_label)[1L]
        else if (!is.null(event) && length(event)) event_info$event <- as.character(event)[1L]
      } else {
        data_raw$.r4vn_score_y <- as.numeric(response)
        if (.r4vn_score_is_binary(response)) {
          event_fit <- if (is.numeric(response) || is.logical(response)) NULL else event
          bi <- .r4vn_score_binary(response, event_fit, "model outcome")
          data_raw$.r4vn_score_y <- bi$y
          event_info <- bi
          if (!is.null(source_event_label) && length(source_event_label)) event_info$event <- as.character(source_event_label)[1L]
          else if (!is.null(event) && length(event)) event_info$event <- as.character(event)[1L]
        } else event_info <- list(event = NULL, levels = NULL)
      }
    }
    if (refprob_supplied) {
      rv0 <- try(eval(refprob_expr, envir = env), silent = TRUE)
      if (!inherits(rv0, "try-error") && is.character(rv0) && length(rv0) == 1L && rv0 %in% names(data_raw)) {
        refprob_full <- as.numeric(data_raw[[rv0]])
        refprob_name <- rv0
      } else if (!inherits(rv0, "try-error") && is.numeric(rv0) && length(rv0) == nrow(data_raw)) {
        refprob_full <- as.numeric(rv0)
      } else {
        stop("With an already fitted model, refprob must be a numeric vector aligned to the model frame or the name of a variable present in that model frame.", call. = FALSE)
      }
      if (any(is.finite(refprob_full) & (refprob_full < 0 | refprob_full > 1)))
        stop("refprob must contain probabilities between 0 and 1 (NA is allowed).", call. = FALSE)
    }
    original_row_index <- seq_len(nrow(data_raw))
    selected <- predictors0
    original_model <- source_model
    select <- "none"
  } else {
    data_input <- if (is.null(data)) .r4vn_score_active_data() else as.data.frame(data)
    outcome_name <- .r4vn_score_resolve_name(out_expr, env)
    pred_expr <- substitute(predictors)
    predictor_spec <- .r4vn_score_predictor_spec(pred_expr, env)
    predictors0 <- predictor_spec$variable
    keep_spec <- !predictor_spec$variable %in% exclude_names
    predictor_spec <- predictor_spec[keep_spec, , drop = FALSE]
    predictors0 <- predictor_spec$variable
    if (!length(predictors0)) stop("No predictors were supplied.", call. = FALSE)
    miss <- setdiff(c(outcome_name, predictors0), names(data_input))
    if (length(miss)) stop("Variables not found in data: ", paste(miss, collapse = ", "), call. = FALSE)
    labels <- setNames(vapply(predictors0, function(v)
      .r4vn_score_label(data_input[[v]], v), character(1L)), predictors0)
    # Honour R4VN vars() declarations: no-prefix/b#. variables are categorical,
    # c./q./f. variables are continuous, and b#. also controls the model reference.
    data_input <- .r4vn_score_prepare_declared(data_input, predictor_spec)
    data_raw_model <- data_input

    y0 <- data_input[[outcome_name]]
    if (time_supplied) time_name <- .r4vn_score_resolve_name(time_expr, env)
    if (!is.null(time_name) && !time_name %in% names(data_input))
      stop("time variable '", time_name, "' was not found in data.", call. = FALSE)
    if (family == "auto") {
      if (!is.null(time_name)) family <- "cox"
      else if (.r4vn_score_is_binary(y0)) family <- "logistic"
      else if (is.numeric(y0) && all(y0[!is.na(y0)] >= 0) &&
               all(abs(y0[!is.na(y0)] - round(y0[!is.na(y0)])) < 1e-8)) family <- "poisson"
      else stop("family='auto' could not determine a supported model. Specify family= explicitly.", call. = FALSE)
    }
    if (family == "cox" && is.null(time_name))
      stop("Cox scorecards require time=followup_variable.", call. = FALSE)

    work <- data_input
    if (family == "logistic") {
      bi <- .r4vn_score_binary(y0, event, outcome_name)
      work$.r4vn_score_y <- bi$y
      event_info <- bi
    } else if (family == "cox") {
      bi <- .r4vn_score_binary(y0, event, outcome_name)
      work$.r4vn_score_y <- bi$y
      work$.r4vn_score_time <- as.numeric(work[[time_name]])
      event_info <- bi
    } else {
      if (.r4vn_score_is_binary(y0)) {
        bi <- .r4vn_score_binary(y0, event, outcome_name)
        work$.r4vn_score_y <- bi$y
        event_info <- bi
      } else {
        work$.r4vn_score_y <- as.numeric(y0)
        event_info <- list(event = NULL, levels = NULL)
      }
    }

    needed <- c(predictors0, ".r4vn_score_y", if (family == "cox") ".r4vn_score_time")
    cc <- stats::complete.cases(work[, needed, drop = FALSE])
    original_row_index <- which(cc)
    work <- work[cc, , drop = FALSE]
    if (nrow(work) < 20L)
      warning("Fewer than 20 complete observations are available for model development.", call. = FALSE)

    if (refprob_supplied) {
      rv0 <- try(eval(refprob_expr, envir = env), silent = TRUE)
      if (!inherits(rv0, "try-error") && is.character(rv0) && length(rv0) == 1L && rv0 %in% names(data_input)) {
        refprob_full <- as.numeric(data_input[[rv0]])
        refprob_name <- rv0
      } else {
        rv <- try(eval(refprob_expr, envir = data_input, enclos = env), silent = TRUE)
        if (inherits(rv, "try-error") || !is.numeric(rv) || length(rv) != nrow(data_input)) {
          stop("refprob must be a numeric vector with one value per original row or a probability-variable name.", call. = FALSE)
        }
        refprob_full <- as.numeric(rv)
        refprob_name <- if (is.symbol(refprob_expr) && as.character(refprob_expr) %in% names(data_input)) as.character(refprob_expr) else NULL
      }
      refprob_full <- refprob_full[cc]
      if (any(is.finite(refprob_full) & (refprob_full < 0 | refprob_full > 1)))
        stop("refprob must contain probabilities between 0 and 1 (NA is allowed).", call. = FALSE)
      if (is.null(refprob_name)) {
        work$.r4vn_score_refprob <- refprob_full
        refprob_name <- ".r4vn_score_refprob"
      }
    }

    force0 <- intersect(force_names, predictors0)
    selected <- .r4vn_score_select(work, family, predictors0, select, force0,
                                    entry, stay, confound, seed)
    if (!length(selected)) stop("No predictors remained after model selection.", call. = FALSE)
    original_model <- .r4vn_score_fit(work, family, selected)
    data_raw <- work
  }

  original_predictor_spec <- .r4vn_score_capture_model_input(data_raw, selected)

  # Build a main-effects score-ready representation. Any continuous
  # categorization therefore remains explicitly testable against the original model.
  tf <- .r4vn_score_transform_fit(data_raw, selected, cuts = cuts,
                                  continuous = continuous, bins = bins)
  score_data <- tf$data
  score_model <- .r4vn_score_fit(score_data, family, selected)
  eff <- .r4vn_score_effect_dictionary(
    score_model, score_data, selected, labels[selected], family,
    riskonly = riskonly, scoreref = scoreref
  )
  if (!isTRUE(simplify) && is.character(points) && length(points) == 1L && points == "auto") points <- "pdo"
  point_build <- .r4vn_score_build_points(
    eff$table, score_data, family, score_model,
    pdo = pdo, points = points, maxscore = maxscore, tolerance = tolerance,
    riskonly = riskonly
  )
  dict <- point_build$dictionary
  model_score <- .r4vn_score_apply_dictionary(score_data, dict, "model_point")
  clinical_score <- .r4vn_score_apply_dictionary(score_data, dict, "clinical_point")
  if (length(unique(clinical_score[is.finite(clinical_score)])) < 2L)
    stop("The clinical point system produces fewer than two distinct total scores. Increase score resolution or revise manual points.", call. = FALSE)
  lp_ready_check <- .r4vn_score_predict_lp(score_model)
  orient <- suppressWarnings(stats::cor(clinical_score, lp_ready_check, use = "complete.obs", method = "spearman"))
  if (isTRUE(riskonly)) {
    if (any(dict$effect < -1e-10, na.rm = TRUE) || any(dict$scoring_ratio < 1 - 1e-10, na.rm = TRUE))
      stop("`riskonly=TRUE` requires every risk-oriented score effect ratio to be >= 1.", call. = FALSE)
    if (any(dict$clinical_point < 0, na.rm = TRUE))
      stop("`riskonly=TRUE` requires a non-negative add-only clinical point system.", call. = FALSE)
  }
  if (is.finite(orient) && orient < 0)
    stop("The supplied clinical points are inversely oriented: higher total points correspond to lower modeled risk. Reverse/revise the manual points so higher score means higher risk before using cutoff-based classification.", call. = FALSE)
  score_range <- .r4vn_score_theoretical_range(dict, "clinical_point")
  attainable_scores <- .r4vn_score_attainable(dict, "clinical_point")
  observed_range <- range(clinical_score, na.rm = TRUE)
  binary_outcome <- .r4vn_score_is_binary(score_data$.r4vn_score_y)

  risk_obj <- NULL
  risk_raw <- NULL
  clinical_model <- NULL
  pred_original <- NULL
  pred_score_model <- NULL
  pred_clinical <- NULL

  if (family == "logistic") {
    pred_original <- as.numeric(stats::predict(original_model, newdata = data_raw, type = "response"))
    pred_score_model <- as.numeric(stats::predict(score_model, newdata = score_data, type = "response"))
    risk_obj <- .r4vn_score_risk_table_binary(score_data$.r4vn_score_y, clinical_score, score_range, riskcut, attainable_scores)
    risk_raw <- .r4vn_score_add_observed_binary(risk_obj$table, score_data$.r4vn_score_y, clinical_score)
    risk_obj$table <- risk_raw
    clinical_model <- risk_obj$fit
    pred_clinical <- as.numeric(stats::predict(clinical_model,
                                               newdata = data.frame(.score = clinical_score),
                                               type = "response"))
  } else if (family == "poisson") {
    pred_original <- as.numeric(stats::predict(original_model, newdata = data_raw, type = "response"))
    pred_score_model <- as.numeric(stats::predict(score_model, newdata = score_data, type = "response"))
    if (binary_outcome) {
      risk_obj <- .r4vn_score_risk_table_binary(score_data$.r4vn_score_y, clinical_score, score_range, riskcut, attainable_scores)
      risk_raw <- .r4vn_score_add_observed_binary(risk_obj$table, score_data$.r4vn_score_y, clinical_score)
      risk_obj$table <- risk_raw
      clinical_model <- risk_obj$fit
      pred_clinical <- as.numeric(stats::predict(clinical_model,
                                                 newdata = data.frame(.score = clinical_score),
                                                 type = "response"))
      pred_original <- pmin(pmax(pred_original, 0), 1)
      pred_score_model <- pmin(pmax(pred_score_model, 0), 1)
    } else {
      risk_obj <- .r4vn_score_risk_table_poisson(score_data$.r4vn_score_y, clinical_score, score_range, attainable_scores)
      risk_raw <- .r4vn_score_add_observed_poisson(risk_obj$table, score_data$.r4vn_score_y, clinical_score)
      risk_obj$table <- risk_raw
      clinical_model <- risk_obj$fit
      pred_clinical <- as.numeric(stats::predict(clinical_model,
                                                 newdata = data.frame(score = clinical_score),
                                                 type = "response"))
    }
  } else if (family == "cox") {
    risk_obj <- .r4vn_score_risk_table_cox(score_data$.r4vn_score_time,
                                           score_data$.r4vn_score_y,
                                           clinical_score, score_range, times, riskcut, attainable_scores)
    risk_raw <- risk_obj$table
    clinical_model <- risk_obj$fit
    times <- risk_obj$times
    cutoff_time <- .r4vn_score_null(cutoff_time, max(times, na.rm = TRUE))
    pred_original <- as.numeric(stats::predict(original_model, newdata = data_raw, type = "lp"))
    pred_score_model <- as.numeric(stats::predict(score_model, newdata = score_data, type = "lp"))
    pred_clinical <- clinical_score
  }

  # Quantify information lost before integer rounding (for example by turning
  # a continuous predictor into bedside categories). This is distinct from the
  # point-rounding loss reported by the automatic point engine.
  representation_loss <- NA_real_
  representation_warning <- NULL
  mo <- .r4vn_score_perf_metric(family, score_data, pred_original)
  ms <- .r4vn_score_perf_metric(family, score_data, pred_score_model)
  if (is.finite(mo) && is.finite(ms)) {
    representation_loss <- mo - ms
    if (representation_loss > tolerance) {
      representation_warning <- paste0(
        "The score-ready predictor representation reduced the primary discrimination metric by ",
        formatC(representation_loss, format = "f", digits = 3),
        ", exceeding tolerance=", formatC(tolerance, format = "f", digits = 3),
        ". Consider clinically specified cuts, more categories, or retaining a more detailed score."
      )
    }
  }

  # Cutoff candidates and the selected final classification rule.
  cutoff_tab <- NULL
  cutoff_perf <- NULL
  selected_cutoff <- NULL
  selected_cutoff_method <- NULL
  refprob_aligned <- refprob_full

  if (cutoff != "none") {
    if (family == "cox") {
      if (cutoff %in% c("prevalence", "refprob")) {
        warning("cutoff='", cutoff,
                "' is not defined for the Cox IPCW workflow; using cutoff='iu'.",
                call. = FALSE)
        cutoff <- "iu"
      }
      cutoff_tab <- .r4vn_score_cutoff_table_cox(
        score_data$.r4vn_score_time, score_data$.r4vn_score_y,
        clinical_score, cutoff_time, risk_raw, riskcut,
        sens, spec, cost_fp, cost_fn, cutoff_value
      )
    } else if (binary_outcome) {
      cutoff_tab <- .r4vn_score_cutoff_table_binary(
        score_data$.r4vn_score_y, clinical_score,
        risk_raw, riskcut, sens, spec, cost_fp, cost_fn,
        cutoff_value, refprob_aligned, refcut
      )
    }

    if (!is.null(cutoff_tab) && nrow(cutoff_tab)) {
      pattern <- switch(
        cutoff,
        iu = "IU",
        youden = "Youden",
        risk = "Clinical risk probability",
        prevalence = "Prevalence probability",
        sens = "Sensitivity",
        spec = "Specificity",
        cost = "cost",
        manual = "Manual",
        refprob = "Reference probability",
        ""
      )
      jj <- grep(pattern, cutoff_tab$Method, ignore.case = TRUE)
      if (!length(jj)) {
        warning("Requested cutoff method could not be calculated; IU/first available cutoff was used.",
                call. = FALSE)
        jj <- grep("IU", cutoff_tab$Method, ignore.case = TRUE)
        if (!length(jj)) jj <- 1L
      }
      sr <- cutoff_tab[jj[[1L]], , drop = FALSE]
      selected_cutoff <- as.numeric(sr$Score_cutoff)
      selected_cutoff_method <- as.character(sr$Method)
      cutoff_tab$Selected <- FALSE
      cutoff_tab$Selected[jj[[1L]]] <- TRUE

      if (family == "cox") {
        cutoff_perf <- .r4vn_score_cut_metrics_cox(
          score_data$.r4vn_score_time, score_data$.r4vn_score_y,
          clinical_score, selected_cutoff, cutoff_time
        )
      } else {
        cutoff_perf <- .r4vn_score_cut_metrics(
          score_data$.r4vn_score_y,
          clinical_score >= selected_cutoff
        )
      }

      if (is.null(riskcut) && !is.null(risk_raw) && nrow(risk_raw)) {
        risk_raw$Risk_group <- ifelse(risk_raw$Score >= selected_cutoff, "High", "Low")
      }
    }
  }

  # Compare original model, score-ready model/model score and clinical score.
  comparison <- NULL
  if (compare) {
    if (family == "logistic" || (family == "poisson" && binary_outcome)) {
      comparison <- .r4vn_score_comparison_binary(
        score_data$.r4vn_score_y,
        pred_original, pred_score_model, pred_clinical, seed
      )
    } else if (family == "cox") {
      comparison <- .r4vn_score_comparison_cox(
        score_data$.r4vn_score_time, score_data$.r4vn_score_y,
        original_model, score_model, clinical_model,
        data_raw, score_data, clinical_score, times
      )
    } else {
      po <- .r4vn_score_poisson_perf(score_data$.r4vn_score_y, pred_original)
      pm <- .r4vn_score_poisson_perf(score_data$.r4vn_score_y, pred_score_model)
      ps <- .r4vn_score_poisson_perf(score_data$.r4vn_score_y, pred_clinical)
      comparison <- data.frame(
        Model = c("Original model", "Model score", "Clinical score"),
        RMSE = c(po$RMSE, pm$RMSE, ps$RMSE),
        MAE = c(po$MAE, pm$MAE, ps$MAE),
        Mean_prediction = c(po$Mean_prediction, pm$Mean_prediction, ps$Mean_prediction),
        stringsAsFactors = FALSE
      )
    }
  }

  time_brier_tab <- if (!is.null(comparison) && family == "cox") attr(comparison, "time_brier") else NULL

  refprob_tab <- NULL
  if (!is.null(refprob_aligned) && !is.null(pred_clinical) && family != "cox") {
    refprob_tab <- .r4vn_score_refprob_table(refprob_aligned, pred_clinical)
  }

  # Publication-oriented and technical model tables.
  model_raw <- .r4vn_score_model_table(original_model, family, source_vcov)
  model_pub <- .r4vn_score_model_publication(
    original_model, family, data_raw, selected, labels[selected], source_vcov
  )

  scorecard_pub <- .r4vn_score_scorecard_display(dict, score_range)
  risk_orientation_pub <- .r4vn_score_orientation_display(dict, eff$effect_measure)
  scorecard_tech <- dict[, c(
    "predictor", "predictor_label", "category",
    "model_effect", "model_ratio", "model_reference",
    "effect", "scoring_ratio", "scoring_reference", "zero_point_category",
    "protective_vs_model_reference", "scoring_rule",
    "model_point", "clinical_point"
  )]

  risk_pub <- NULL
  if (risktable && !is.null(risk_raw)) {
    if (family == "poisson" && !binary_outcome) {
      risk_pub <- data.frame(
        Score = risk_raw$Score,
        `Expected value (95% CI)` = paste0(
          .r4vn_score_num(risk_raw$Predicted_mean, 3), " (",
          .r4vn_score_num(risk_raw$CI_low, 3), "\u2013",
          .r4vn_score_num(risk_raw$CI_high, 3), ")"
        ),
        check.names = FALSE
      )
    } else {
      risk_pub <- data.frame(
        Score = risk_raw$Score,
        `Predicted risk (95% CI)` = paste0(
          .r4vn_score_pct(risk_raw$Predicted_risk, 1), " (",
          .r4vn_score_pct(risk_raw$CI_low, 1), "\u2013",
          .r4vn_score_pct(risk_raw$CI_high, 1), ")"
        ),
        check.names = FALSE
      )
      if ("Time" %in% names(risk_raw))
        risk_pub <- cbind(Time = risk_raw$Time, risk_pub)
      if ("Risk_group" %in% names(risk_raw))
        risk_pub$`Risk group` <- risk_raw$Risk_group
    }
  }

  cutoff_pub <- cutoff_tab
  if (!is.null(cutoff_pub) && nrow(cutoff_pub)) {
    for (nm in intersect(c("Risk_threshold", "Sensitivity", "Specificity",
                           "Youden", "IU", "AUC"), names(cutoff_pub))) {
      cutoff_pub[[nm]] <- round(cutoff_pub[[nm]], 3)
    }
  }

  cutoff_perf_pub <- cutoff_perf
  if (!is.null(cutoff_perf_pub) && nrow(cutoff_perf_pub)) {
    cutoff_perf_pub$`Estimate (95% CI)` <- ifelse(
      is.finite(cutoff_perf_pub$CI_low),
      paste0(
        .r4vn_score_num(cutoff_perf_pub$Estimate, 3), " (",
        .r4vn_score_num(cutoff_perf_pub$CI_low, 3), "\u2013",
        .r4vn_score_num(cutoff_perf_pub$CI_high, 3), ")"
      ),
      .r4vn_score_num(cutoff_perf_pub$Estimate, 3)
    )
    keepn <- c("Measure", "Estimate (95% CI)",
               if ("Time" %in% names(cutoff_perf_pub)) "Time")
    cutoff_perf_pub <- cutoff_perf_pub[, keepn, drop = FALSE]
  }

  decision_data <- NULL
  calibration_data <- NULL
  roc_data <- NULL
  if (family == "logistic" || (family == "poisson" && binary_outcome)) {
    if (decision)
      decision_data <- .r4vn_score_decision_curve(
        score_data$.r4vn_score_y, pred_original, pred_clinical
      )
    if (calibration)
      calibration_data <- .r4vn_score_calibration_data(
        score_data$.r4vn_score_y, pred_clinical
      )
    roc_data <- .r4vn_score_roc_table(score_data$.r4vn_score_y, clinical_score)
  } else if (family == "cox" && !is.null(cutoff_time)) {
    roc_data <- .r4vn_score_time_roc(
      score_data$.r4vn_score_time, score_data$.r4vn_score_y,
      clinical_score, cutoff_time
    )
  }

  rebuild_args <- NULL
  if (!supplied_model) {
    rebuild_args <- list(
      outcome = outcome_name,
      predictors = if (!is.null(predictor_spec)) .r4vn_score_rebuild_predictor_tokens(predictor_spec) else predictors0,
      data = NULL,
      family = family,
      event = event_info$event,
      time = time_name,
      times = times,
      cutoff_time = cutoff_time,
      select = select,
      force = force_names,
      exclude = NULL,
      entry = entry,
      stay = stay,
      confound = confound,
      cuts = cuts,
      continuous = continuous,
      bins = bins,
      points = points,
      pdo = pdo,
      maxscore = maxscore,
      simplify = simplify,
      tolerance = tolerance,
      riskonly = riskonly,
      scoreref = scoreref,
      cutoff = cutoff,
      riskcut = riskcut,
      cutoff_value = cutoff_value,
      sens = sens,
      spec = spec,
      cost_fp = cost_fp,
      cost_fn = cost_fn,
      refprob = refprob_name,
      refcut = refcut,
      validate = "none",
      bootstrap = 0L,
      validation = NULL,
      risktable = risktable,
      compare = compare,
      calibration = calibration,
      decision = FALSE,
      plot = FALSE,
      show = FALSE,
      console = FALSE,
      seed = seed,
      ai = FALSE
    )
  }

  ans <- list(
    call = match.call(),
    family = family,
    n = nrow(score_data),
    binary_outcome = binary_outcome,
    event = event_info$event,
    outcome = outcome_name,
    time = time_name,
    times = times,
    cutoff_time = cutoff_time,
    candidate_predictors = predictors0,
    selected_predictors = selected,
    predictor_spec = predictor_spec,
    original_predictor_spec = original_predictor_spec,
    transform = tf$spec,
    labels = labels[selected],
    score_method = point_build$method,
    effect_measure = eff$effect_measure,
    riskonly = isTRUE(riskonly),
    scoreref = scoreref,
    pdo = pdo,
    scale_B = point_build$scale_B,
    simplification_warning = point_build$warning,
    simplification_loss = point_build$performance_loss,
    representation_warning = representation_warning,
    representation_loss = representation_loss,
    score_range = score_range,
    attainable_scores = attainable_scores,
    observed_score_range = observed_range,
    selected_cutoff = selected_cutoff,
    selected_cutoff_method = selected_cutoff_method,
    riskcut = riskcut,
    models = list(
      original = original_model,
      score_model = score_model,
      clinical = clinical_model
    ),
    scores = list(
      model = model_score,
      clinical = clinical_score
    ),
    predictions = list(
      original = pred_original,
      model_score = pred_score_model,
      clinical = pred_clinical
    ),
    tables = list(
      model = model_pub,
      model_raw = model_raw,
      scorecard = scorecard_pub,
      risk_orientation = risk_orientation_pub,
      scorecard_technical = scorecard_tech,
      risk = risk_pub,
      risk_raw = risk_raw,
      cutoff = cutoff_pub,
      cutoff_raw = cutoff_tab,
      cutoff_performance = cutoff_perf_pub,
      cutoff_performance_raw = cutoff_perf,
      comparison = comparison,
      time_brier = time_brier_tab,
      reference_probability = refprob_tab,
      validation = NULL,
      external_validation = NULL
    ),
    plot_enabled = isTRUE(plot),
    plot_data = list(
      risk = risk_raw,
      roc = roc_data,
      calibration = calibration_data,
      decision = decision_data,
      distribution = data.frame(
        Score = clinical_score,
        Outcome = score_data$.r4vn_score_y
      )
    ),
    .dictionary = dict,
    .development_data = score_data,
    .development_data_raw = data_raw,
    .raw_original_data = if (!is.null(data_raw_model)) data_raw_model else data_raw,
    .original_row_index = original_row_index,
    .refprob = refprob_aligned,
    .rebuild_args = rebuild_args
  )
  class(ans) <- c("r4vn_tabscore", "list")
  ans$plots <- .r4vn_score_available_plots(ans)
  ans$plot_titles <- setNames(
    vapply(ans$plots, function(z) .r4vn_score_plot_title(ans, z), character(1)),
    ans$plots
  )

  if (!is.null(validation)) {
    ev <- try(.r4vn_score_external_validation(ans, validation), silent = TRUE)
    if (inherits(ev, "try-error")) {
      warning("External validation could not be completed: ", as.character(ev),
              call. = FALSE)
    } else ans$tables$external_validation <- ev
  }

  if (validate == "bootstrap") {
    if (is.null(rebuild_args)) {
      warning(
        "Full-pipeline bootstrap validation is unavailable when tabscore() is called with an already fitted model. Re-run tabscore() from outcome + candidate predictors to validate development.",
        call. = FALSE
      )
    } else {
      ans$tables$validation <- .r4vn_score_bootstrap_validation(
        ans, as.integer(bootstrap), seed
      )
    }
  }

  # Flat, non-NULL publication tables are kept separately so they can be
  # passed directly to R4VN::tabexport(), which accepts lists of data frames.
  pub <- list(
    `Final model` = ans$tables$model,
    `Risk-oriented score coding` = ans$tables$risk_orientation,
    `Clinical scorecard` = ans$tables$scorecard,
    `Score to risk` = ans$tables$risk,
    `Cutoff selection` = ans$tables$cutoff,
    `Cutoff performance` = ans$tables$cutoff_performance,
    `Model versus score` = ans$tables$comparison,
    `Time-specific Brier` = ans$tables$time_brier,
    `Reference probability` = ans$tables$reference_probability,
    `Internal validation` = ans$tables$validation,
    `External validation` = ans$tables$external_validation
  )
  ans$publication_tables <- Filter(function(z) is.data.frame(z) || is.matrix(z), pub)

  if (!identical(ai, FALSE)) {
    if (exists("aiask", mode = "function", inherits = TRUE)) {
      aa <- try({
        if (isTRUE(ai))
          get("aiask", mode = "function", inherits = TRUE)(ans)
        else
          get("aiask", mode = "function", inherits = TRUE)(ans, api = ai)
      }, silent = TRUE)
      if (!inherits(aa, "try-error")) ans$ai <- aa else ans$ai <- NULL
    } else {
      warning("ai= was requested but R4VN aiask() is not available.", call. = FALSE)
    }
  }

  if (isTRUE(plot) && interactive() && length(ans$plots)) {
    try(graphics::plot(ans, which = "all"), silent = TRUE)
  }

  shown <- FALSE
  if (isTRUE(show)) {
    hp <- try(.r4vn_score_show_html(ans, show = TRUE), silent = TRUE)
    if (!inherits(hp, "try-error")) {
      ans$file <- hp
      shown <- TRUE
    }
  }
  if (isTRUE(console) || (!shown && isTRUE(show))) print(ans)
  invisible(ans)
}

#' Predict scores and risk from an R4VN scorecard
#'
#' @param object A `r4vn_tabscore` object.
#' @param newdata New data containing all final scorecard predictors.
#' @param type Output type: `all`, `score`, `risk`, `group`, `model`, `model_lp`,
#'   or `model_score`. `model` returns original-model response prediction for
#'   logistic/Poisson and linear predictor for Cox. `model_lp` always returns the
#'   original-model linear predictor.
#' @param times Cox prediction horizon(s). Defaults to the horizons stored in the
#'   scorecard.
#' @param ... Reserved.
#' @return A numeric vector, matrix, or data frame depending on `type`.
#' @method predict r4vn_tabscore
#' @export
predict.r4vn_tabscore <- function(object, newdata, type = c("all", "score", "risk", "group",
                                                            "model", "model_lp", "model_score"),
                                  times = NULL, ...) {
  type <- match.arg(type)
  newdata <- as.data.frame(newdata)
  sd <- .r4vn_score_transform_apply(newdata, object$transform)
  score <- .r4vn_score_apply_dictionary(sd, object$.dictionary, "clinical_point")
  mscore <- .r4vn_score_apply_dictionary(sd, object$.dictionary, "model_point")

  if (type == "score") return(score)
  if (type == "model_score") return(mscore)
  if (type == "model_lp") {
    tp <- if (object$family == "cox") "lp" else "link"
    nd0 <- .r4vn_score_apply_model_input(newdata, object$original_predictor_spec)
    return(as.numeric(stats::predict(object$models$original, newdata = nd0, type = tp)))
  }
  if (type == "model") {
    tp <- if (object$family == "cox") "lp" else "response"
    nd0 <- .r4vn_score_apply_model_input(newdata, object$original_predictor_spec)
    z <- as.numeric(stats::predict(object$models$original, newdata = nd0, type = tp))
    if (object$family == "poisson" && object$binary_outcome) z <- pmin(pmax(z, 0), 1)
    return(z)
  }

  if (object$family == "logistic" || (object$family == "poisson" && object$binary_outcome)) {
    risk <- as.numeric(stats::predict(object$models$clinical,
                                     newdata = data.frame(.score = score),
                                     type = "response"))
    group <- if (!is.null(object$riskcut)) {
      .r4vn_score_risk_group(risk, object$riskcut)
    } else if (!is.null(object$selected_cutoff)) {
      ifelse(score >= object$selected_cutoff, "High", "Low")
    } else rep(NA_character_, length(score))
    if (type == "risk") return(risk)
    if (type == "group") return(group)
    return(data.frame(Score = score, Predicted_risk = risk, Risk_group = group,
                      stringsAsFactors = FALSE))
  }

  if (object$family == "poisson") {
    risk <- as.numeric(stats::predict(object$models$clinical,
                                     newdata = data.frame(score = score),
                                     type = "response"))
    if (type == "risk") return(risk)
    if (type == "group") return(rep(NA_character_, length(score)))
    return(data.frame(Score = score, Predicted_value = risk))
  }

  tt <- .r4vn_score_null(times, object$times)
  if (is.null(tt) || !length(tt))
    stop("Cox risk prediction requires times= or stored prediction horizons.", call. = FALSE)
  dtmp <- data.frame(.r4vn_score_time = 1, .r4vn_score_y = 1, .score = score)
  risk <- .r4vn_score_cox_risk(object$models$clinical,
                               data.frame(.score = score), tt)
  colnames(risk) <- paste0("Risk_", tt)
  if (type == "risk") return(risk)
  if (type == "group") {
    if (is.null(object$riskcut)) return(matrix(NA_character_, nrow(risk), ncol(risk),
                                               dimnames = dimnames(risk)))
    g <- apply(risk, 2L, .r4vn_score_risk_group, cuts = object$riskcut)
    return(g)
  }
  out <- data.frame(Score = score)
  cbind(out, as.data.frame(risk, check.names = FALSE))
}

#' @export
print.r4vn_tabscore <- function(x, ...) {
  cat("R4VN Scorecard\n")
  cat("Model: ", x$family, "\n", sep = "")
  cat("N: ", x$n, "\n", sep = "")
  cat("Predictors: ", paste(x$selected_predictors, collapse = ", "), "\n", sep = "")
  cat("Clinical score range: ", x$score_range[["min"]], "\u2013", x$score_range[["max"]], "\n", sep = "")
  if (isTRUE(x$riskonly)) cat("Score orientation: risk-only, add-only (all clinical points >= 0; score effects >= 1).\n")
  if (!is.null(x$selected_cutoff))
    cat("Selected cutoff: \u2265", x$selected_cutoff, " (", x$selected_cutoff_method, ")\n", sep = "")
  if (!is.null(x$representation_loss) && is.finite(x$representation_loss))
    cat("Discrimination loss from score-ready categorization: ", format(round(x$representation_loss, 4), nsmall = 4), "\n", sep = "")
  if (!is.null(x$simplification_loss) && is.finite(x$simplification_loss))
    cat("Additional discrimination loss from integer points: ", format(round(x$simplification_loss, 4), nsmall = 4), "\n", sep = "")
  if (!is.null(x$representation_warning)) cat("Note: ", x$representation_warning, "\n", sep = "")
  if (!is.null(x$simplification_warning)) cat("Note: ", x$simplification_warning, "\n", sep = "")
  invisible(x)
}

#' @export
summary.r4vn_tabscore <- function(object, ...) {
  list(
    family = object$family,
    n = object$n,
    predictors = object$selected_predictors,
    score_range = object$score_range,
    riskonly = object$riskonly,
    effect_measure = object$effect_measure,
    representation_loss = object$representation_loss,
    point_simplification_loss = object$simplification_loss,
    cutoff = object$selected_cutoff,
    cutoff_method = object$selected_cutoff_method,
    comparison = object$tables$comparison,
    validation = object$tables$validation
  )
}

.r4vn_score_external_validation <- function(object, validation) {
  d <- as.data.frame(validation)
  if (is.null(object$outcome) || !object$outcome %in% names(d)) {
    stop("External validation needs the original outcome variable '",
         .r4vn_score_null(object$outcome, "<unknown>"), "'.", call. = FALSE)
  }
  if (object$family == "cox") {
    if (is.null(object$time) || !object$time %in% names(d))
      stop("External Cox validation needs time variable '", object$time, "'.", call. = FALSE)
    bi <- .r4vn_score_binary(d[[object$outcome]], object$event, object$outcome)
    ok <- stats::complete.cases(d[, unique(c(object$selected_predictors,
                                              object$outcome, object$time)), drop = FALSE])
    dv <- d[ok, , drop = FALSE]
    ev <- bi$y[ok]
    tm <- as.numeric(dv[[object$time]])
    lp <- stats::predict(object, dv, type = "model_lp")
    sc <- stats::predict(object, dv, type = "score")
    c1 <- .r4vn_score_cindex(tm, ev, lp)
    c2 <- .r4vn_score_cindex(tm, ev, sc)
    hz <- object$cutoff_time
    b1 <- b2 <- NA_real_
    if (!is.null(hz) && is.finite(hz)) {
      nd0 <- .r4vn_score_apply_model_input(dv, object$original_predictor_spec)
      pr1 <- try(.r4vn_score_cox_risk(object$models$original, nd0, hz)[, 1L], silent = TRUE)
      pr2 <- try(stats::predict(object, dv, type = "risk", times = hz)[, 1L], silent = TRUE)
      if (!inherits(pr1, "try-error")) b1 <- .r4vn_score_ipcw_brier(tm, ev, pr1, hz)
      if (!inherits(pr2, "try-error")) b2 <- .r4vn_score_ipcw_brier(tm, ev, pr2, hz)
    }
    out <- data.frame(
      Model = c("Original model", "Clinical score"),
      C_index = c(c1[1], c2[1]),
      CI_low = c(c1[2], c2[2]),
      CI_high = c(c1[3], c2[3]),
      IPCW_Brier = c(b1, b2),
      Time = if (is.null(hz)) NA_real_ else hz,
      Sensitivity = NA_real_, Specificity = NA_real_, PPV = NA_real_, NPV = NA_real_, Accuracy = NA_real_,
      N = nrow(dv), stringsAsFactors = FALSE
    )
    if (!is.null(object$selected_cutoff) && is.finite(object$selected_cutoff) &&
        !is.null(hz) && is.finite(hz)) {
      cm <- try(.r4vn_score_cut_metrics_cox(tm, ev, sc, object$selected_cutoff, hz), silent = TRUE)
      if (!inherits(cm, "try-error")) {
        getm_cox_external <- function(prefix) {
          jj <- grep(paste0("^", prefix), cm$Measure, ignore.case = TRUE)
          if (!length(jj)) return(NA_real_)
          as.numeric(cm$Estimate[jj[[1L]]])
        }
        out$Sensitivity[2L] <- getm_cox_external("Sensitivity")
        out$Specificity[2L] <- getm_cox_external("Specificity")
        out$PPV[2L] <- getm_cox_external("PPV")
        out$NPV[2L] <- getm_cox_external("NPV")
        out$Accuracy[2L] <- getm_cox_external("Accuracy")
      }
    }
    return(out)
  }

  y0 <- d[[object$outcome]]
  if (object$binary_outcome) {
    bi <- .r4vn_score_binary(y0, object$event, object$outcome)
    ok <- stats::complete.cases(d[, unique(c(object$selected_predictors,
                                              object$outcome)), drop = FALSE])
    dv <- d[ok, , drop = FALSE]
    y <- bi$y[ok]
    po <- stats::predict(object, dv, type = "model")
    ps <- stats::predict(object, dv, type = "risk")
    sc <- stats::predict(object, dv, type = "score")
    a1 <- .r4vn_score_auc_ci(y, po)
    a2 <- .r4vn_score_auc_ci(y, ps)
    ca <- .r4vn_score_calibration(y, po)
    cs <- .r4vn_score_calibration(y, ps)
    out <- data.frame(
      Model = c("Original model", "Clinical score"),
      AUC = c(a1[1], a2[1]),
      AUC_low = c(a1[2], a2[2]),
      AUC_high = c(a1[3], a2[3]),
      Brier = c(mean((y - po)^2, na.rm = TRUE), mean((y - ps)^2, na.rm = TRUE)),
      Calibration_intercept = c(ca[["intercept"]], cs[["intercept"]]),
      Calibration_slope = c(ca[["slope"]], cs[["slope"]]),
      Sensitivity = NA_real_, Specificity = NA_real_, PPV = NA_real_, NPV = NA_real_, Accuracy = NA_real_,
      N = nrow(dv), stringsAsFactors = FALSE
    )
    if (!is.null(object$selected_cutoff) && is.finite(object$selected_cutoff)) {
      cm <- .r4vn_score_cut_metrics(y, sc >= object$selected_cutoff)
      getm_binary_external <- function(nm) cm$Estimate[match(nm, cm$Measure)]
      out$Sensitivity[2L] <- getm_binary_external("Sensitivity")
      out$Specificity[2L] <- getm_binary_external("Specificity")
      out$PPV[2L] <- getm_binary_external("PPV")
      out$NPV[2L] <- getm_binary_external("NPV")
      out$Accuracy[2L] <- getm_binary_external("Accuracy")
    }
    return(out)
  }

  ok <- stats::complete.cases(d[, unique(c(object$selected_predictors,
                                            object$outcome)), drop = FALSE])
  dv <- d[ok, , drop = FALSE]
  y <- as.numeric(dv[[object$outcome]])
  po <- stats::predict(object, dv, type = "model")
  ps <- stats::predict(object, dv, type = "risk")
  p1 <- .r4vn_score_poisson_perf(y, po)
  p2 <- .r4vn_score_poisson_perf(y, ps)
  data.frame(Model = c("Original model", "Clinical score"),
             RMSE = c(p1$RMSE, p2$RMSE), MAE = c(p1$MAE, p2$MAE),
             Mean_prediction = c(p1$Mean_prediction, p2$Mean_prediction),
             N = nrow(dv), stringsAsFactors = FALSE)
}

#' Plot an R4VN scorecard
#'
#' @description
#' Draw one or all publication-oriented scorecard graphics using base R only.
#' With `which="all"` (the default), every available graph is drawn in sequence;
#' in RStudio the back/forward arrows in the Plots pane can be used to review the
#' complete plot history. The same available graphics are embedded automatically
#' in the HTML Viewer when the original `tabscore()` call used `plot=TRUE`.
#'
#' @param x A `r4vn_tabscore` object.
#' @param which Plot type: `all`, `risk`, `roc`, `calibration`, `decision`, or
#'   `distribution`. `all` draws every plot that is available for the fitted
#'   model family.
#' @param title Optional custom title when one plot is requested. With
#'   `which="all"`, each plot keeps its own descriptive title.
#' @param font_family Base-R graphics font family. Default `"sans"` is used to
#'   keep Viewer, browser and RStudio rendering consistent without another
#'   graphics dependency.
#' @param ... Additional arguments are reserved.
#' @return For one plot, invisibly returns its plotted data. With `which="all"`,
#'   invisibly returns a named list containing the data for every graph drawn.
#' @examples
#' set.seed(23)
#' d <- data.frame(
#'   age = rnorm(160, 50, 11),
#'   smoke = factor(rbinom(160, 1, .30), 0:1, c("No", "Yes"))
#' )
#' d$event <- rbinom(160, 1, plogis(-3.5 + .045*d$age + .7*(d$smoke == "Yes")))
#' z <- tabscore(event, c(age, smoke), data=d,
#'               validate="none", plot=TRUE, show=FALSE)
#' plot(z, which="risk")
#' plot(z, which="roc")
#' \donttest{
#' plot(z) # all available plots, one after another
#' }
#' @export
plot.r4vn_tabscore <- function(
    x,
    which = c("all", "risk", "roc", "calibration", "decision", "distribution"),
    title = NULL,
    font_family = "sans",
    ...
) {
  if (!inherits(x, "r4vn_tabscore"))
    stop("`x` must be created by tabscore().", call. = FALSE)
  which <- match.arg(which)
  available <- .r4vn_score_available_plots(x)
  if (!length(available)) stop("No plot data are available for this scorecard.", call. = FALSE)

  if (which == "all") {
    out <- lapply(available, function(nm) {
      .r4vn_score_plot_draw(x, nm, font_family = font_family)
    })
    names(out) <- available
    return(invisible(out))
  }

  if (!which %in% available)
    stop("No plot data are available for '", which, "'.", call. = FALSE)
  .r4vn_score_plot_draw(x, which, title = title, font_family = font_family)
}

Try the R4VN package in your browser

Any scripts or data that you put into this service are public.

R4VN documentation built on Sept. 30, 2026, 5:13 p.m.