Nothing
# 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("&", "&", x, fixed = TRUE)
x <- gsub("<", "<", x, fixed = TRUE)
x <- gsub(">", ">", x, fixed = TRUE)
x <- gsub('"', """, 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> | n = ", x$n,
" | 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)
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.