Nothing
R/tabsurv.R
In R4VN: Health Data Analysis and Publication-Ready Reporting
Defines functions
print.r4vn_surv
.r4vn_tabsurv_plot
.r4vn_tabsurv_contract
.r4vn_tabsurv_hierarchical_result
.r4vn_tabsurv_profile
.r4vn_tabsurv_by_info
.r4vn_tabsurv_core
tabsurv
Documented in
tabsurv
#' Comprehensive Survival Analysis Table
#'
#' Performs descriptive survival analysis, Kaplan-Meier/Aalen-Johansen estimates,
#' optional life tables, cumulative incidence at selected times, incidence rate,
#' log-rank tests, Cox regression, proportional-hazards diagnostics, RMST,
#' competing-risk Fine-Gray models, and counting-process/recurrent-event Cox models.
#'
#' @param time Follow-up or stop-time variable, supplied without quotes.
#' @param event Event/status variable, supplied without quotes.
#' @param vars Optional predictor specification created by `vars()`.
#' @param by Optional grouping variable for survival curves and comparisons.
#' Hierarchical syntax is supported: in `by = vars(province, sex, treatment)`,
#' treatment is the innermost curve/comparison group and province > sex are
#' ordered outer strata.
#' @param data Optional data frame. When omitted, active R4VN data are used.
#' @param failure Value of `event` representing the event of interest. For a
#' binary event it defaults to the second factor level or larger numeric value.
#' @param compete Optional competing-event value(s). When supplied, `risk = TRUE`
#' uses the Aalen-Johansen cumulative incidence function.
#' @param id Optional subject identifier for counting-process/recurrent data.
#' @param start Optional start/entry time. When supplied, `time` is treated as stop time.
#' @param unit Optional display unit such as "day", "month", or "year".
#' @param followup Estimate median follow-up using reverse Kaplan-Meier when possible.
#' With `report = "auto"`, this is enabled unless explicitly set to `FALSE`.
#' @param km Fit Kaplan-Meier (ordinary survival) or Aalen-Johansen (competing risks).
#' @param lifetable Show a detailed life table at every observed time. The
#' default is `FALSE`. For ordinary survival, the table reports numbers at
#' risk, events, censoring, conditional survival, cumulative Kaplan-Meier
#' survival, cumulative risk, standard error, and confidence limits. With
#' competing risks, it reports the corresponding Aalen-Johansen event-history
#' table and cumulative incidence.
#' @param at Optional time points for survival/risk/rate summaries. With
#' `report = "auto"` or `"full"`, three representative follow-up times are
#' selected automatically when `at` is omitted.
#' @param risk Report cumulative risk at `at`. For ordinary survival this is 1-S(t);
#' with competing risks it is the cumulative incidence function.
#' @param cuminc Optional numeric time points at which cumulative incidence is
#' required, for example `cuminc = c(6, 12, 24)`. This directly activates
#' cumulative-risk output without also requiring `risk = TRUE`. Ordinary
#' survival uses 1-KM; competing-risk analysis uses the Aalen-Johansen CIF.
#' @param rate `FALSE`, `TRUE`, `"overall"`, `"cumulative"`, `"interval"`,
#' or `"all"`. `"all"` reports overall, cumulative, and interval-specific
#' rates. The automatic profile uses the overall incidence rate.
#' @param scale Rate multiplier, e.g. 100 for events per 100 person-time units.
#' @param logrank Perform a log-rank test when `by` is supplied and no competing risk exists.
#' @param rr,rd Compare cumulative risks between two `by` groups using approximate
#' risk ratio or risk difference inference based on survival-estimate standard errors.
#' @param irr Compare incidence rates between two `by` groups.
#' @param cox Fit crude Cox models for variables in `vars`.
#' @param adjusted FALSE/NULL, TRUE (adjust each focal predictor for all other focal
#' predictors), or a `vars()`/character set of adjustment covariates.
#' @param multi FALSE/NULL, TRUE (all `vars` in one model), or a `vars()`/character
#' set defining the final multivariable Cox model.
#' @param strata Optional stratification variable for Cox regression.
#' @param cluster Optional clustering variable for robust Cox variance.
#' @param frailty Optional shared-frailty variable. Do not combine with `cluster`.
#' @param finegray Fit a Fine-Gray subdistribution hazards model when `compete` is supplied.
#' @param recurrent FALSE/TRUE or "ag". TRUE is Andersen-Gill and requires
#' `id` and `start`. Automatic profiles suppress ordinary KM/RMST modules for
#' recurrent-event data unless the user explicitly requests them.
#' @param rmst Compute restricted mean survival time.
#' @param tau Restriction time for RMST. Defaults to the largest common curve time.
#' @param ph Test the proportional-hazards assumption with `cox.zph()` for the final Cox model.
#' @param interaction Optional `vars(a, b)` containing exactly two variables to include
#' their interaction in the final multivariable Cox model.
#' @param superby Optional outer subgroup variable retained for backward compatibility.
#' For new code, multiple ordered outer strata can be supplied directly in
#' `by = vars(stratum1, stratum2, group)`.
#' @param ci Confidence level, default 0.95.
#' @param digit,p_digit,effect_digit Display digits.
#' @param missing Show missing/exclusion information when printing.
#' @param plot Draw a survival/CIF curve using `gsurv()` after analysis. In the
#' automatic profiles, the plot includes confidence limits, the log-rank
#' p-value when available, and a number-at-risk table.
#' @param title Optional title.
#' @param show Logical; open the formatted result in the Viewer. Default `TRUE`.
#' @param console Logical; also print the traditional result in the Console. Default `FALSE`.
#' @param ai Prepare a compact de-identified interpretation payload in `$ai_text`.
#' @param ties Cox tie method: "efron", "breslow", or "exact".
#' @param report Reporting profile: `"auto"` (context-sensitive comprehensive
#' output), `"brief"` (descriptive survival summary), `"full"` (all valid
#' modules), or `"custom"` (backward-compatible concise defaults plus explicitly requested
#' modules).
#' @param plot_args Named list of additional arguments passed to `gsurv()`.
#' @param interpretation Add a cautious, deterministic interpretation table.
#' The default is `FALSE`; use `interpretation = TRUE` when narrative output
#' is wanted.
#' @param export Optional export format accepted by `tabexport()`, such as
#' `"docx"`, `"xlsx"`, or `"html"`.
#' @param file Optional export filename. Its extension may also determine the
#' export format.
#' @param open Open the exported file when supported.
#' @param strict If `TRUE`, an unavailable optional module stops the analysis.
#' The default `FALSE` keeps the main report and records a warning instead.
#'
#' @return An object of class `r4vn_surv`. Backward-compatible components are
#' retained, with a consistent reporting contract in `$descriptive`,
#' `$estimates`, `$tests`, `$diagnostics`, `$interpretation`, `$tables`,
#' `$plots`, `$models`, `$metadata`, and `$call`.
#' @export
#'
#' @examples
#' if (requireNamespace("survival", quietly = TRUE)) {
#' # Reproducible two-group data from the survival package.
#' d <- survival::lung
#' d$death <- as.integer(d$status == 2)
#' d$group <- factor(d$sex, levels = c(1, 2),
#' labels = c("Male", "Female"))
#'
#' # 1. Complete two-group report. This includes the log-rank test.
#' km <- tabsurv(
#' time, death, by = group, data = d, failure = 1,
#' unit = "day", at = c(90, 180, 365, 540),
#' report = "auto", plot = FALSE, show = FALSE
#' )
#' km$logrank
#' km$tests$logrank
#' km$logrank$p
#'
#' \donttest{
#' # 2. Cumulative incidence at 6, 12, and 24 months.
#' d$month <- d$time / 30.4375
#' ci_month <- tabsurv(
#' month, death, by = group, data = d, failure = 1,
#' cuminc = c(6, 12, 24), report = "custom", show = FALSE
#' )
#' ci_month$cuminc
#'
#' # 3. Detailed life table and interpretation are both opt-in.
#' km_detail <- tabsurv(
#' time, death, by = group, data = d, failure = 1,
#' at = c(90, 180, 365, 540),
#' lifetable = TRUE, interpretation = TRUE, show = FALSE
#' )
#' head(km_detail$lifetable)
#' km_detail$interpretation
#'
#' # 4. A compact KM plus log-rank analysis without automatic extras.
#' km_simple <- tabsurv(
#' time, death, by = group, data = d, failure = 1,
#' report = "custom", km = TRUE, logrank = TRUE,
#' plot = FALSE, show = FALSE
#' )
#'
#' # 5. Explicit two-group effect measures and RMST.
#' km_compare <- tabsurv(
#' time, death, by = group, data = d, failure = 1,
#' at = c(90, 180, 365, 540), risk = TRUE,
#' rr = TRUE, rd = TRUE, rate = "all", irr = TRUE,
#' rmst = TRUE, tau = 365, show = FALSE
#' )
#' km_compare$risk_compare
#' km_compare$irr
#' km_compare$rmst
#'
#' # 6. Publication graphs, including risk tables, are documented in ?gsurv.
#' # Keeping graphics out of this example also keeps tabsurv() examples fast
#' # and executable on non-interactive CRAN check devices.
#'
#' # 7. With competing risks, use the Aalen-Johansen CIF, not 1-KM.
#' set.seed(2026)
#' n <- 180
#' t1 <- rexp(n, 0.07)
#' t2 <- rexp(n, 0.05)
#' tc <- runif(n, 4, 30)
#' tm <- pmin(t1, t2, tc)
#' dcr <- data.frame(
#' time = tm,
#' status = ifelse(tm == t1, 1L, ifelse(tm == t2, 2L, 0L)),
#' group = factor(rep(c("A", "B"), each = n / 2))
#' )
#' cif <- tabsurv(
#' time, status, by = group, data = dcr,
#' failure = 1, compete = 2, cuminc = c(6, 12, 24),
#' report = "custom", show = FALSE
#' )
#' cif$cuminc
#' # See ?gsurv for publication CIF graphs and risk tables.
#' }
#' }
tabsurv <- function(time, event, vars = NULL, by = NULL, data = NULL,
failure = NULL, compete = NULL,
id = NULL, start = NULL, unit = NULL,
followup = NULL, km = NULL, lifetable = FALSE, at = NULL,
risk = NULL, cuminc = NULL, rate = NULL, scale = 100,
logrank = NULL, rr = NULL, rd = NULL, irr = NULL,
cox = NULL, adjusted = NULL, multi = NULL,
strata = NULL, cluster = NULL, frailty = NULL,
finegray = NULL, recurrent = FALSE,
rmst = NULL, tau = NULL,
ph = NULL, interaction = FALSE, superby = NULL,
ci = 0.95, digit = 2, p_digit = 3, effect_digit = 2,
missing = FALSE, plot = NULL,
title = NULL, show = TRUE, console = FALSE, ai = FALSE,
ties = c("efron", "breslow", "exact"),
report = c("auto", "brief", "full", "custom"),
plot_args = list(), interpretation = FALSE,
export = NULL, file = NULL, open = FALSE,
strict = FALSE) {
.r4vn_surv_require()
call <- match.call()
env <- parent.frame()
d <- .r4vn_surv_data(data)
report <- match.arg(report)
ties <- match.arg(ties)
if (!is.list(plot_args) ||
(length(plot_args) && (is.null(names(plot_args)) || any(!nzchar(names(plot_args)))))) {
stop("`plot_args` must be a named list.", call. = FALSE)
}
if (!is.logical(interpretation) || length(interpretation) != 1L || is.na(interpretation)) {
stop("`interpretation` must be TRUE or FALSE.", call. = FALSE)
}
if (!is.logical(lifetable) || length(lifetable) != 1L || is.na(lifetable)) {
stop("`lifetable` must be TRUE or FALSE.", call. = FALSE)
}
if (!is.null(cuminc) && (!is.numeric(cuminc) || !length(cuminc) ||
any(!is.finite(cuminc)) || any(cuminc < 0))) {
stop("`cuminc` must contain non-negative finite time points.", call. = FALSE)
}
if (!is.null(cuminc)) cuminc <- sort(unique(cuminc))
if (!is.logical(strict) || length(strict) != 1L || is.na(strict)) {
stop("`strict` must be TRUE or FALSE.", call. = FALSE)
}
by_expr <- substitute(by)
by_info <- .r4vn_tabsurv_by_info(by_expr, d, env)
time_name <- .r4vn_surv_name(substitute(time), d, "time")
provided <- c(
followup = !missing(followup), km = !missing(km), at = !missing(at),
risk = !missing(risk), rate = !missing(rate), logrank = !missing(logrank),
rr = !missing(rr), rd = !missing(rd), irr = !missing(irr),
cox = !missing(cox), adjusted = !missing(adjusted), multi = !missing(multi),
finegray = !missing(finegray), rmst = !missing(rmst), ph = !missing(ph),
plot = !missing(plot)
)
profile <- .r4vn_tabsurv_profile(
report = report, provided = provided,
values = list(
followup = followup, km = km, at = at, risk = risk, rate = rate,
logrank = logrank, rr = rr, rd = rd, irr = irr, cox = cox,
adjusted = adjusted, multi = multi, finegray = finegray, rmst = rmst,
ph = ph, plot = plot, vars = vars
),
data = d, time_name = time_name, by_name = by_info$by,
competing = !is.null(compete), recurrent = recurrent
)
for (nm in names(profile)) assign(nm, profile[[nm]])
if (!is.null(cuminc)) risk <- TRUE
core_call <- call
core_call[[1L]] <- quote(.r4vn_tabsurv_core)
core_call$data <- quote(.r4vn_tabsurv_data)
core_call$followup <- followup
core_call$km <- km
core_call$lifetable <- lifetable
core_call$at <- at
core_call$risk <- risk
core_call$cuminc <- cuminc
core_call$rate <- rate
core_call$logrank <- logrank
core_call$rr <- rr
core_call$rd <- rd
core_call$irr <- irr
core_call$cox <- cox
if (!isTRUE(provided[["adjusted"]]) || is.logical(adjusted)) core_call$adjusted <- adjusted
if (!isTRUE(provided[["multi"]]) || is.logical(multi)) core_call$multi <- multi
core_call$finegray <- finegray
core_call$rmst <- rmst
core_call$ph <- ph
core_call$plot <- plot
core_call$show <- FALSE
core_call$console <- FALSE
core_call$ties <- ties
core_call$report <- report
core_call$strict <- strict
core_call$plot_args <- NULL
core_call$interpretation <- NULL
core_call$export <- NULL
core_call$file <- NULL
core_call$open <- NULL
if (!is.null(by_info$spec)) {
spec <- by_info$spec
outer <- spec$strata
super_expr <- substitute(superby)
if (!.r4vn_expr_is_null(super_expr)) {
explicit <- .r4vn_resolve_name_spec(
super_expr, d, env, "superby", allow_null = TRUE, multiple = TRUE
)
outer <- unique(c(explicit, outer))
}
outer <- setdiff(outer, spec$by)
if (length(outer)) {
ids <- .r4vn_strata_indices(d, outer)
if (!length(ids)) stop("No complete strata are available for survival analysis.", call. = FALSE)
labels <- vapply(ids, function(idx) .r4vn_stratum_label(d, outer, idx), character(1))
core_call$by <- as.name(spec$by)
core_call$superby <- NULL
core_call$plot <- FALSE
core_call$ai <- FALSE
results <- lapply(ids, function(idx) {
ee <- new.env(parent = env)
ee$.r4vn_tabsurv_core <- .r4vn_tabsurv_core
ee$.r4vn_tabsurv_data <- d[idx, , drop = FALSE]
z <- eval(core_call, envir = ee)
z$call <- call
.r4vn_tabsurv_contract(z, interpretation = interpretation)
})
names(results) <- labels
out <- .r4vn_tabsurv_hierarchical_result(
results, labels, spec, outer, d, call, title, missing, ci, scale, report
)
if (isTRUE(ai)) {
out$ai_text <- lapply(results, function(z) {
tryCatch(.r4vn_surv_ai_text(z, effect_digit, p_digit),
error = function(e) NULL)
})
}
} else {
core_call$by <- as.name(spec$by)
core_call$superby <- NULL
ee <- new.env(parent = env)
ee$.r4vn_tabsurv_core <- .r4vn_tabsurv_core
ee$.r4vn_tabsurv_data <- d
out <- eval(core_call, envir = ee)
}
} else {
ee <- new.env(parent = env)
ee$.r4vn_tabsurv_core <- .r4vn_tabsurv_core
ee$.r4vn_tabsurv_data <- d
out <- eval(core_call, envir = ee)
}
out$call <- call
out <- .r4vn_tabsurv_contract(out, interpretation = interpretation)
out <- .r4vn_tabsurv_plot(
out, plot = plot, plot_args = plot_args, report = report, show = show
)
if (!is.null(export) || !is.null(file)) {
out$export <- survexport(
out, export = export, file = file, open = open,
title = title %||% "R4VN survival analysis"
)
}
.r4vn_show(out, show = show, console = console)
}
.r4vn_tabsurv_core <- function(time, event, vars = NULL, by = NULL, data = NULL,
failure = NULL, compete = NULL,
id = NULL, start = NULL, unit = NULL,
followup = TRUE, km = TRUE, lifetable = FALSE, at = NULL,
risk = FALSE, cuminc = NULL, rate = FALSE, scale = 100,
logrank = TRUE, rr = FALSE, rd = FALSE, irr = FALSE,
cox = FALSE, adjusted = FALSE, multi = FALSE,
strata = NULL, cluster = NULL, frailty = NULL,
finegray = FALSE, recurrent = FALSE,
rmst = FALSE, tau = NULL,
ph = FALSE, interaction = FALSE, superby = NULL,
ci = 0.95, digit = 2, p_digit = 3, effect_digit = 2,
missing = FALSE, plot = FALSE,
title = NULL, show = FALSE, console = FALSE, ai = FALSE,
ties = c("efron", "breslow", "exact"),
report = "custom", strict = FALSE) {
.r4vn_surv_require()
call <- match.call()
env <- parent.frame()
data <- .r4vn_surv_data(data)
ties <- match.arg(ties)
if (!is.numeric(ci) || length(ci) != 1L || is.na(ci) || ci <= 0 || ci >= 1) stop("`ci` must be between 0 and 1.", call. = FALSE)
if (!is.numeric(scale) || length(scale) != 1L || is.na(scale) || scale <= 0) stop("`scale` must be a positive number.", call. = FALSE)
if (!is.null(at)) {
if (!is.numeric(at) || any(!is.finite(at)) || any(at < 0)) stop("`at` must contain non-negative finite time points.", call. = FALSE)
at <- sort(unique(at))
}
if (!is.null(cuminc)) {
if (!is.numeric(cuminc) || !length(cuminc) || any(!is.finite(cuminc)) ||
any(cuminc < 0)) {
stop("`cuminc` must contain non-negative finite time points.", call. = FALSE)
}
cuminc <- sort(unique(cuminc))
risk <- TRUE
}
time_name <- .r4vn_surv_name(substitute(time), data, "time")
event_name <- .r4vn_surv_name(substitute(event), data, "event")
by_name <- .r4vn_surv_name(substitute(by), data, "by", TRUE)
id_name <- .r4vn_surv_name(substitute(id), data, "id", TRUE)
start_name <- .r4vn_surv_name(substitute(start), data, "start", TRUE)
strata_name0 <- .r4vn_surv_name(substitute(strata), data, "strata", TRUE)
cluster_name0 <- .r4vn_surv_name(substitute(cluster), data, "cluster", TRUE)
frailty_name0 <- .r4vn_surv_name(substitute(frailty), data, "frailty", TRUE)
superby_name <- .r4vn_surv_name(substitute(superby), data, "superby", TRUE)
if (!is.null(cluster_name0) && !is.null(frailty_name0)) stop("Use either `cluster` or `frailty`, not both in the same Cox model.", call. = FALSE)
if (!isFALSE(recurrent)) {
rec <- if (isTRUE(recurrent)) "ag" else tolower(as.character(recurrent)[1L])
if (!identical(rec, "ag")) stop("This implementation currently supports recurrent = TRUE or recurrent = 'ag' (Andersen-Gill).", call. = FALSE)
if (is.null(id_name) || is.null(start_name)) stop("Andersen-Gill recurrent-event analysis requires both `id` and `start`.", call. = FALSE)
if (is.null(cluster_name0)) cluster_name0 <- id_name
cox <- TRUE
multi <- if (isFALSE(multi)) TRUE else multi
}
status <- .r4vn_surv_event_status(data[[event_name]], failure, compete)
.r4vn_surv_validate_time(data[[time_name]], status$event, if (is.null(start_name)) NULL else data[[start_name]])
focal_spec <- .r4vn_surv_as_spec(vars, data, "vars", TRUE)
adjusted_spec <- NULL
if (!isFALSE(adjusted) && !is.null(adjusted)) {
adjusted_spec <- if (isTRUE(adjusted)) TRUE else .r4vn_surv_as_spec(adjusted, data, "adjusted", FALSE)
}
multi_spec <- NULL
if (!isFALSE(multi) && !is.null(multi)) {
multi_spec <- if (isTRUE(multi)) TRUE else .r4vn_surv_as_spec(multi, data, "multi", FALSE)
}
interaction_vars <- NULL
if (!isFALSE(interaction) && !is.null(interaction)) {
isp <- .r4vn_surv_as_spec(interaction, data, "interaction", FALSE)
if (nrow(isp) != 2L) stop("`interaction` must specify exactly two variables, for example `vars(treatment, sex)`.", call. = FALSE)
interaction_vars <- isp$variable
multi <- if (isFALSE(multi)) TRUE else multi
multi_spec <- if (is.null(multi_spec)) TRUE else multi_spec
}
all_spec <- .r4vn_surv_merge_specs(
focal_spec,
if (isTRUE(adjusted_spec)) NULL else adjusted_spec,
if (isTRUE(multi_spec)) NULL else multi_spec
)
if (isTRUE(adjusted_spec) || isTRUE(multi_spec)) all_spec <- .r4vn_surv_merge_specs(all_spec, focal_spec)
internal <- .r4vn_surv_internal_data(data, all_spec)
d <- internal$data
map <- internal$map
d$.time <- data[[time_name]]
d$.event <- status$event
d$.status_ms <- status$status_ms
if (!is.null(start_name)) d$.start <- data[[start_name]]
if (!is.null(by_name)) d$.by <- if (is.factor(data[[by_name]])) droplevels(data[[by_name]]) else factor(data[[by_name]], levels = unique(data[[by_name]][!is.na(data[[by_name]])]))
if (!is.null(id_name)) d$.id <- data[[id_name]]
if (!is.null(strata_name0)) d$.strata <- data[[strata_name0]]
if (!is.null(cluster_name0)) d$.cluster <- data[[cluster_name0]]
if (!is.null(frailty_name0)) d$.frailty <- data[[frailty_name0]]
if (!is.null(superby_name)) d$.superby <- if (is.factor(data[[superby_name]])) droplevels(data[[superby_name]]) else factor(data[[superby_name]])
# Use separate analysis sets. Descriptive survival estimates must not lose
# participants merely because a Cox covariate is missing.
competing <- !is.null(compete)
if (competing && !is.null(start_name) && is.null(id_name)) {
stop("Competing-risk start-stop data require `id`.", call. = FALSE)
}
desc_required <- c(".time", ".event")
if (!is.null(start_name)) desc_required <- c(desc_required, ".start")
if (!is.null(by_name) && (km || isTRUE(lifetable) || risk || !is.null(cuminc) ||
!identical(rate, FALSE) || logrank || rr || rd || irr || rmst ||
isTRUE(plot))) {
desc_required <- c(desc_required, ".by")
}
if (!is.null(start_name) && !is.null(id_name)) desc_required <- c(desc_required, ".id")
cc_desc <- .r4vn_surv_complete(d, desc_required)
d_desc <- cc_desc$data
if (!nrow(d_desc)) stop("No complete observations remain for survival-time analysis.", call. = FALSE)
if (".by" %in% names(d_desc)) d_desc$.by <- droplevels(d_desc$.by)
# Cox/Fine-Gray models use their own base set. Predictor-specific missingness
# is left to coxph()/finegray(), so crude models are not forced to use the
# complete-case sample of every other predictor.
model_required <- c(".time", ".event")
if (!is.null(start_name)) model_required <- c(model_required, ".start")
if (!is.null(strata_name0)) model_required <- c(model_required, ".strata")
if (!is.null(cluster_name0)) model_required <- c(model_required, ".cluster")
if (!is.null(frailty_name0)) model_required <- c(model_required, ".frailty")
if (!is.null(superby_name)) model_required <- c(model_required, ".superby")
if ((!isFALSE(recurrent) || (isTRUE(finegray) && !is.null(start_name))) && !is.null(id_name)) model_required <- c(model_required, ".id")
cc_model <- .r4vn_surv_complete(d, model_required)
d_model <- cc_model$data
if (!nrow(d_model)) stop("No complete observations remain for the requested model structure.", call. = FALSE)
if (!is.null(map) && nrow(map)) {
for (i in seq_len(nrow(map))) if (map$type[i] == "categorical") d_model[[map$internal[i]]] <- droplevels(d_model[[map$internal[i]]])
}
if (".superby" %in% names(d_model)) d_model$.superby <- droplevels(d_model$.superby)
overview <- .r4vn_surv_overview(d_desc)
descriptive <- .r4vn_surv_group_overview(d_desc)
followup_table <- if (isTRUE(followup)) .r4vn_surv_reverse_followup(d_desc, ci) else NULL
fit <- curve <- median_table <- at_table <- risk_table <- life_table <- NULL
logrank_table <- NULL
risk_times <- if (!is.null(cuminc)) cuminc else at
if (km || risk || rmst || isTRUE(lifetable) || !is.null(at) ||
!is.null(cuminc) || isTRUE(plot)) {
rhs <- if (".by" %in% names(d_desc)) ".by" else "1"
if (!competing) {
sf <- .r4vn_surv_formula(".start" %in% names(d_desc), rhs)
fit <- survival::survfit(sf, data = d_desc, conf.int = ci,
id = if (".start" %in% names(d_desc) && ".id" %in% names(d_desc)) d_desc$.id else NULL)
curve <- .r4vn_surv_curve_km(fit)
median_table <- .r4vn_surv_median(fit)
if (!is.null(at)) {
at_table <- .r4vn_surv_km_at(fit, at)
}
if (isTRUE(risk) && !is.null(risk_times)) {
rs <- if (!is.null(at_table) && identical(risk_times, at)) {
at_table
} else .r4vn_surv_km_at(fit, risk_times)
if (nrow(rs)) {
rs$risk <- 1 - rs$surv
rs$risk_lower <- 1 - rs$upper
rs$risk_upper <- 1 - rs$lower
rs$se <- rs$std.err
risk_table <- rs[, intersect(c("group", "time", "n.risk", "risk", "se", "risk_lower", "risk_upper", "surv", "lower", "upper"), names(rs)), drop = FALSE]
}
}
} else {
sf <- .r4vn_surv_ms_formula(".start" %in% names(d_desc), rhs)
fit <- survival::survfit(sf, data = d_desc, conf.int = ci, id = if (".id" %in% names(d_desc)) d_desc$.id else NULL)
curve <- .r4vn_surv_curve_cif(fit, "failure")
if (isTRUE(risk) && !is.null(risk_times)) {
risk_table <- .r4vn_surv_cif_at(curve, d_desc, risk_times)
names(risk_table)[names(risk_table) == "cif"] <- "risk"
# Use the Aalen-Johansen standard error; keep a CI-based fallback for
# older survival objects that do not expose std.err.
zcrit <- stats::qnorm(1 - (1 - ci) / 2)
risk_table$se <- risk_table$std.err
bad_se <- !is.finite(risk_table$se)
risk_table$se[bad_se] <- (risk_table$upper[bad_se] - risk_table$lower[bad_se]) / (2 * zcrit)
risk_table$risk_lower <- risk_table$lower
risk_table$risk_upper <- risk_table$upper
}
}
if (isTRUE(lifetable) && !is.null(curve) && nrow(curve)) {
life_table <- .r4vn_surv_lifetable(curve, competing = competing)
}
}
if (isTRUE(logrank) && ".by" %in% names(d_desc) && ".start" %in% names(d_desc) && !competing) {
warning("Log-rank testing is not reported for start-stop/counting-process data; use the Cox model for comparison.", call. = FALSE)
}
if (isTRUE(logrank) && ".by" %in% names(d_desc) && nlevels(d_desc$.by) > 1L && !competing && !(".start" %in% names(d_desc))) {
lf <- .r4vn_surv_formula(FALSE, ".by")
lr <- survival::survdiff(lf, data = d_desc)
df <- max(1L, length(lr$n) - 1L)
logrank_table <- data.frame(chisq = unname(lr$chisq), df = df,
p = stats::pchisq(lr$chisq, df, lower.tail = FALSE), stringsAsFactors = FALSE)
}
rate_table <- irr_table <- NULL
if (!identical(rate, FALSE)) {
mode <- if (is.character(rate)) tolower(rate[1L]) else TRUE
allowed_rate <- c("overall", "cumulative", "interval", "all")
if (is.character(mode) && !mode %in% allowed_rate) {
stop("`rate` must be FALSE, TRUE, 'overall', 'cumulative', 'interval', or 'all'.", call. = FALSE)
}
rate_table <- .r4vn_surv_rate_table(d_desc, at, mode, scale, ci)
if (isTRUE(irr)) {
if (!".by" %in% names(d_desc)) stop("`irr = TRUE` requires `by`.", call. = FALSE)
irr_table <- tryCatch(
.r4vn_surv_compare_rates(rate_table, ci),
error = function(e) {
if (isTRUE(strict)) stop(e)
warning("Incidence-rate comparison was not available: ", conditionMessage(e), call. = FALSE)
NULL
}
)
}
}
risk_compare <- NULL
if (isTRUE(rr) || isTRUE(rd)) {
if (is.null(risk_times) || !length(risk_times)) {
stop("`rr`/`rd` require `at` or `cuminc` time points.", call. = FALSE)
}
if (!".by" %in% names(d_desc)) stop("`rr`/`rd` require `by`.", call. = FALSE)
if (is.null(risk_table)) stop("Risk estimates are unavailable.", call. = FALSE)
risk_compare <- tryCatch(
.r4vn_surv_compare_risk(risk_table, ci, rr, rd),
error = function(e) {
if (isTRUE(strict)) stop(e)
warning("Cumulative-risk comparison was not available: ", conditionMessage(e), call. = FALSE)
NULL
}
)
}
rmst_result <- NULL
if (isTRUE(rmst)) {
if (competing) stop("RMST in this function is for ordinary survival. For competing risks use CIF/Fine-Gray results.", call. = FALSE)
if (is.null(fit)) {
rhs <- if (".by" %in% names(d_desc)) ".by" else "1"
fit <- survival::survfit(.r4vn_surv_formula(".start" %in% names(d_desc), rhs), data = d_desc, conf.int = ci)
}
if (is.null(tau)) {
if (is.null(fit$strata)) tau <- max(fit$time, na.rm = TRUE) else {
ends <- cumsum(as.integer(fit$strata)); starts <- c(1L, head(ends, -1L) + 1L)
tau <- min(vapply(Map(seq.int, starts, ends), function(ii) max(fit$time[ii], na.rm = TRUE), numeric(1)))
}
}
if (!is.numeric(tau) || length(tau) != 1L || !is.finite(tau) || tau <= 0) stop("`tau` must be one positive finite time point.", call. = FALSE)
rmst_result <- tryCatch(
.r4vn_surv_rmst(fit, tau, ci),
error = function(e) {
if (isTRUE(strict)) stop(e)
warning("RMST was not available: ", conditionMessage(e), call. = FALSE)
NULL
}
)
}
# Cox models -------------------------------------------------------------
cox_result <- list(crude = NULL, adjusted = NULL, multi = NULL,
interaction = NULL, diagnostics = NULL,
crude_fits = list(), adjusted_fits = list(), multi_fit = NULL)
need_cox <- isTRUE(cox) || !isFALSE(adjusted) || !isFALSE(multi) || isTRUE(ph) || !isFALSE(recurrent)
if (need_cox) {
if (is.null(focal_spec) || !nrow(focal_spec)) stop("Cox analysis requires `vars`.", call. = FALSE)
focal_map <- map[match(focal_spec$variable, map$variable), , drop = FALSE]
adj_obj <- if (isTRUE(adjusted_spec)) TRUE else if (is.null(adjusted_spec)) NULL else adjusted_spec
mul_obj <- if (isTRUE(multi_spec) || (isTRUE(multi) && is.null(multi_spec))) TRUE else if (is.null(multi_spec)) NULL else multi_spec
cox_result <- tryCatch(
.r4vn_surv_cox_models(
d_model, focal_map, map, adj_obj, mul_obj,
start = ".start" %in% names(d_model),
strata_name = if (".strata" %in% names(d_model)) ".strata" else NULL,
cluster_name = if (".cluster" %in% names(d_model)) ".cluster" else NULL,
frailty_name = if (".frailty" %in% names(d_model)) ".frailty" else NULL,
ci = ci, ties = ties, interaction_vars = interaction_vars
),
error = function(e) {
if (isTRUE(strict)) stop(e)
warning("Cox analysis was not available: ", conditionMessage(e), call. = FALSE)
list(crude = NULL, adjusted = NULL, multi = NULL,
interaction = NULL, diagnostics = NULL,
crude_fits = list(), adjusted_fits = list(), multi_fit = NULL)
}
)
if (!isTRUE(cox)) cox_result$crude <- NULL
}
ph_table <- NULL
if (isTRUE(ph)) {
ph_table <- tryCatch({
ph_fit <- cox_result$multi_fit
if (is.null(ph_fit)) {
focal_map <- map[match(focal_spec$variable, map$variable), , drop = FALSE]
ph_fit <- .r4vn_surv_fit_cox(
d_model, focal_map, ".start" %in% names(d_model),
if (".strata" %in% names(d_model)) ".strata" else NULL,
if (".cluster" %in% names(d_model)) ".cluster" else NULL,
if (".frailty" %in% names(d_model)) ".frailty" else NULL,
ties = ties
)
if (is.null(cox_result$multi_fit)) cox_result$multi_fit <- ph_fit
}
.r4vn_surv_ph(ph_fit)
}, error = function(e) {
if (isTRUE(strict)) stop(e)
warning("PH test unavailable: ", conditionMessage(e), call. = FALSE)
NULL
})
}
fg_result <- NULL
if (isTRUE(finegray)) {
if (!competing) stop("`finegray = TRUE` requires `compete`.", call. = FALSE)
fg_spec <- if (!is.null(multi_spec) && !isTRUE(multi_spec)) multi_spec else focal_spec
if (is.null(fg_spec)) stop("`finegray = TRUE` requires `vars`.", call. = FALSE)
fg_map <- map[match(fg_spec$variable, map$variable), , drop = FALSE]
fg_result <- tryCatch(
.r4vn_surv_finegray(d_model, fg_map, ".start" %in% names(d_model), ci, ties),
error = function(e) {
if (isTRUE(strict)) stop(e)
warning("Fine-Gray analysis was not available: ", conditionMessage(e), call. = FALSE)
NULL
}
)
}
# Subgroup final Cox models ---------------------------------------------
subgroup <- NULL
if (".superby" %in% names(d_model)) {
if (is.null(focal_spec) || !nrow(focal_spec)) stop("`superby` requires `vars` for subgroup Cox models.", call. = FALSE)
subgroup <- lapply(levels(d_model$.superby), function(g) {
z <- d_model[d_model$.superby == g, , drop = FALSE]
focal_map <- map[match(focal_spec$variable, map$variable), , drop = FALSE]
fit <- try(.r4vn_surv_fit_cox(z, focal_map, ".start" %in% names(z),
if (".strata" %in% names(z)) ".strata" else NULL,
if (".cluster" %in% names(z)) ".cluster" else NULL,
if (".frailty" %in% names(z)) ".frailty" else NULL,
ties = ties), silent = TRUE)
if (inherits(fit, "try-error")) return(list(group = g, fit = NULL, table = NULL))
list(group = g, fit = fit, table = .r4vn_surv_extract_cox(fit, focal_map, ci))
})
names(subgroup) <- levels(d_model$.superby)
}
# Compact analysis data kept for reusable plotting/risk tables.
analysis_data <- d_desc[, intersect(c(".time", ".event", ".start", ".by", ".id", ".status_ms"), names(d_desc)), drop = FALSE]
out <- list(
title = title,
overview = overview,
descriptive = descriptive,
followup = followup_table,
median = median_table,
lifetable = life_table,
at = at_table,
risk = if (isTRUE(risk) || isTRUE(rr) || isTRUE(rd) || !is.null(cuminc)) risk_table else NULL,
cuminc = if (isTRUE(risk) || isTRUE(rr) || isTRUE(rd) || !is.null(cuminc)) risk_table else NULL,
rate = rate_table,
risk_compare = risk_compare,
irr = irr_table,
logrank = logrank_table,
rmst = rmst_result,
cox = cox_result,
ph = ph_table,
finegray = fg_result,
subgroup = subgroup,
fit = fit,
curve = curve,
analysis_data = analysis_data,
metadata = list(time = time_name, event = event_name, by = by_name,
time_label = .r4vn_surv_label(data, time_name),
event_label = .r4vn_surv_label(data, event_name),
by_label = .r4vn_surv_label(data, by_name), id = id_name,
start = start_name, unit = unit, failure = status$failure, compete = compete,
vars = focal_spec, ci = ci, scale = scale, at = at,
cuminc = cuminc, competing = competing,
excluded = cc_desc$excluded, excluded_survival = cc_desc$excluded,
excluded_model_base = cc_model$excluded, missing = missing,
recurrent = recurrent, report = report,
lifetable = isTRUE(lifetable)),
call = call
)
class(out) <- "r4vn_surv"
if (isTRUE(ai)) out$ai_text <- .r4vn_surv_ai_text(out, effect_digit, p_digit)
out
}
.r4vn_tabsurv_by_info <- function(expr, data, env) {
if (.r4vn_expr_is_null(expr)) return(list(by = NULL, spec = NULL))
if (.r4vn_is_vars_spec_expr(expr, data, env)) {
spec <- .r4vn_by_spec(expr, data, env, allow_null = FALSE)
return(list(by = spec$by, spec = spec))
}
list(by = .r4vn_surv_name(expr, data, "by", TRUE), spec = NULL)
}
.r4vn_tabsurv_profile <- function(report, provided, values, data, time_name,
by_name = NULL, competing = FALSE,
recurrent = FALSE) {
n_vars <- if (is.null(values$vars)) 0L else if (inherits(values$vars, "r4vn_vars")) {
nrow(.r4vn_resolve_vars(values$vars, data = data, default_type = "auto", strict = TRUE))
} else if (is.character(values$vars)) length(unique(values$vars)) else 1L
has_vars <- n_vars > 0L
has_by <- !is.null(by_name)
n_groups <- if (has_by) {
length(unique(as.character(data[[by_name]][!is.na(data[[by_name]])])))
} else 0L
two_groups <- n_groups == 2L
legacy <- list(
followup = TRUE, km = TRUE, at = NULL, risk = FALSE, rate = FALSE,
logrank = TRUE, rr = FALSE, rd = FALSE, irr = FALSE,
cox = FALSE, adjusted = FALSE, multi = FALSE, finegray = FALSE,
rmst = FALSE, ph = FALSE, plot = FALSE
)
brief <- utils::modifyList(legacy, list(plot = TRUE))
auto <- list(
followup = TRUE, km = TRUE,
at = .r4vn_surv_auto_times(data[[time_name]]),
risk = TRUE, rate = "overall",
logrank = has_by && !competing,
rr = two_groups, rd = two_groups, irr = two_groups,
cox = has_vars, adjusted = FALSE, multi = n_vars > 1L,
finegray = competing && has_vars,
rmst = two_groups && !competing,
ph = has_vars, plot = TRUE
)
full <- utils::modifyList(auto, list(
adjusted = n_vars > 1L,
multi = has_vars,
rmst = has_by && !competing,
rate = "all"
))
defaults <- switch(report, custom = legacy, brief = brief, auto = auto, full = full)
out <- defaults
for (nm in names(out)) {
if (isTRUE(provided[[nm]]) && !is.null(values[[nm]])) out[[nm]] <- values[[nm]]
}
if (isTRUE(provided[["at"]])) out$at <- values$at
if (!isFALSE(recurrent)) {
out$cox <- TRUE
out$multi <- TRUE
for (nm in c("km", "risk", "logrank", "rr", "rd", "rmst", "plot", "finegray")) {
if (!isTRUE(provided[[nm]])) out[[nm]] <- FALSE
}
if (!isTRUE(provided[["at"]])) out$at <- NULL
}
if (competing) {
if (!isTRUE(provided[["logrank"]])) out$logrank <- FALSE
if (!isTRUE(provided[["rmst"]])) out$rmst <- FALSE
}
out
}
.r4vn_tabsurv_hierarchical_result <- function(results, labels, spec, outer,
data, call, title, missing,
ci, scale, report) {
n_total <- sum(vapply(results, function(z) z$overview$N[1L], numeric(1)), na.rm = TRUE)
ev_total <- sum(vapply(results, function(z) z$overview$Events[1L], numeric(1)), na.rm = TRUE)
cens_total <- sum(vapply(results, function(z) z$overview$Censored[1L], numeric(1)), na.rm = TRUE)
pt_total <- sum(vapply(results, function(z) z$overview$`Person-time`[1L], numeric(1)), na.rm = TRUE)
overview <- data.frame(
N = n_total, Events = ev_total, Censored = cens_total,
`Event percent` = if (n_total > 0L) 100 * ev_total / n_total else NA_real_,
`Person-time` = pt_total, check.names = FALSE
)
out <- list(
title = title, overview = overview,
hierarchical_results = results, strata_results = results,
strata_labels = labels,
hierarchical_by = list(
all = spec$all, strata = outer, by = spec$by,
all_labels = vapply(spec$all, function(nm) .r4vn_variable_label(data, nm), character(1)),
strata_labels = vapply(outer, function(nm) .r4vn_variable_label(data, nm), character(1)),
by_label = .r4vn_variable_label(data, spec$by)
),
metadata = list(
hierarchical = TRUE, by = spec$by, strata = outer, all_by = spec$all,
missing = missing, ci = ci, scale = scale, report = report
),
call = call
)
class(out) <- c("r4vn_surv_hierarchical", "r4vn_surv")
out
}
.r4vn_tabsurv_contract <- function(x, interpretation = FALSE) {
if (!inherits(x, "r4vn_surv")) return(x)
if (!is.null(x$hierarchical_results)) {
x$descriptive <- do.call(rbind, lapply(seq_along(x$hierarchical_results), function(i) {
z <- x$hierarchical_results[[i]]$descriptive
if (is.null(z) || !nrow(z)) return(NULL)
data.frame(Stratum = names(x$hierarchical_results)[i], z,
stringsAsFactors = FALSE, check.names = FALSE)
}))
x$estimates <- lapply(x$hierarchical_results, `[[`, "estimates")
x$tests <- lapply(x$hierarchical_results, `[[`, "tests")
x$diagnostics <- lapply(x$hierarchical_results, `[[`, "diagnostics")
x$models <- lapply(x$hierarchical_results, `[[`, "models")
if (isTRUE(interpretation)) {
x$interpretation <- do.call(rbind, lapply(seq_along(x$hierarchical_results), function(i) {
z <- x$hierarchical_results[[i]]$interpretation
if (is.null(z) || !nrow(z)) return(NULL)
data.frame(Stratum = names(x$hierarchical_results)[i], z,
stringsAsFactors = FALSE, check.names = FALSE)
}))
} else x$interpretation <- NULL
} else {
x$estimates <- list(
survival_at = x$at, life_table = x$lifetable,
cumulative_risk = x$risk, cumulative_incidence = x$cuminc,
incidence_rate = x$rate, rmst = x$rmst,
cox = x$cox[c("crude", "adjusted", "multi", "interaction")],
finegray = if (is.null(x$finegray)) NULL else x$finegray$table
)
x$tests <- list(
logrank = x$logrank, risk_comparison = x$risk_compare,
incidence_rate_ratio = x$irr
)
x$diagnostics <- list(
cox = if (is.null(x$cox)) NULL else x$cox$diagnostics,
proportional_hazards = x$ph
)
x$models <- list(
survival = x$fit,
crude_cox = if (is.null(x$cox)) list() else x$cox$crude_fits,
adjusted_cox = if (is.null(x$cox)) list() else x$cox$adjusted_fits,
multivariable_cox = if (is.null(x$cox)) NULL else x$cox$multi_fit,
finegray = if (is.null(x$finegray)) NULL else x$finegray$fit
)
if (isTRUE(interpretation)) {
x$interpretation <- .r4vn_surv_interpretation(x)
} else x$interpretation <- NULL
}
x$tables <- surv_tables(x)
x
}
.r4vn_tabsurv_plot <- function(x, plot = FALSE, plot_args = list(),
report = "custom", show = TRUE) {
if (!isTRUE(plot)) return(x)
gf <- get0("gsurv", mode = "function", inherits = TRUE)
if (is.null(gf)) {
warning("`plot = TRUE` requested but `gsurv()` is not available.", call. = FALSE)
return(x)
}
defaults <- if (identical(report, "custom")) list() else list(
ci = TRUE, pvalue = TRUE, risk_table = TRUE, median = FALSE
)
args <- utils::modifyList(defaults, plot_args)
if ("x" %in% names(args)) args$x <- NULL
if (is.null(args$show)) args$show <- isTRUE(show)
if (!isTRUE(show)) args$show <- FALSE
draw_one <- function(z, label = NULL) tryCatch({
args2 <- args
if (is.null(args2$legend) && !is.null(z$metadata$by_label)) {
args2$legend <- z$metadata$by_label
}
if (!is.null(label) && !is.null(args2$file)) {
ext <- tools::file_ext(args2$file)
stem <- if (nzchar(ext)) {
substr(args2$file, 1L, nchar(args2$file) - nchar(ext) - 1L)
} else args2$file
suffix <- gsub("[^A-Za-z0-9]+", "-", label)
suffix <- gsub("(^-+|-+$)", "", suffix)
args2$file <- paste0(stem, "-", suffix, if (nzchar(ext)) paste0(".", ext) else "")
}
do.call(gf, c(list(x = z), args2))
},
error = function(e) {
warning("Survival graph was not available: ", conditionMessage(e), call. = FALSE)
NULL
}
)
if (!is.null(x$hierarchical_results)) {
x$graphs <- Map(draw_one, x$hierarchical_results, names(x$hierarchical_results))
names(x$graphs) <- names(x$hierarchical_results)
x$plots <- x$graphs
} else {
x$graph <- draw_one(x)
x$plots <- list(survival = x$graph)
}
x
}
#' @export
print.r4vn_surv <- function(x, ...) {
if (!is.null(x$hierarchical_results)) {
if (!is.null(x$title) && nzchar(x$title)) cat(x$title, "\n", sep = "")
cat("R4VN hierarchical survival analysis\n")
cat(strrep("-", 58), "\n", sep = "")
cat("Hierarchy: ", paste(x$hierarchical_by$strata_labels %||%
x$hierarchical_by$strata, collapse = " > "),
" | Within-stratum group: ", x$hierarchical_by$by_label %||%
x$hierarchical_by$by, "\n", sep = "")
for (i in seq_along(x$hierarchical_results)) {
lab <- names(x$hierarchical_results)[i] %||% paste0("Stratum ", i)
cat("\n=== ", lab, " ===\n", sep = "")
print(x$hierarchical_results[[i]])
}
return(invisible(x))
}
if (!is.null(x$title) && nzchar(x$title)) cat(x$title, "\n", sep = "")
cat("R4VN survival analysis\n")
cat(strrep("-", 58), "\n", sep = "")
ov <- x$overview
cat(sprintf("N: %d | Events: %d | Censored: %d | Person-time: %.2f\n",
ov$N[1L], ov$Events[1L], ov$Censored[1L], ov$`Person-time`[1L]))
if (!is.null(x$descriptive) && nrow(x$descriptive)) {
cat("\nOutcome by ", x$metadata$by_label %||% "group", "\n", sep = "")
print(x$descriptive, row.names = FALSE)
}
if (isTRUE(x$metadata$missing) || x$metadata$excluded_survival > 0L || x$metadata$excluded_model_base > 0L) {
cat(sprintf("Excluded from survival summaries because time/event/group fields were incomplete: %d\n", x$metadata$excluded_survival))
if (x$metadata$excluded_model_base != x$metadata$excluded_survival) {
cat(sprintf("Excluded from the model base because structural model fields were incomplete: %d\n", x$metadata$excluded_model_base))
}
}
if (!is.null(x$followup) && nrow(x$followup)) {
cat("Median follow-up: ", .r4vn_surv_ci_text(x$followup$median[1], x$followup$lower[1], x$followup$upper[1], 2), "\n", sep = "")
}
if (!is.null(x$median) && nrow(x$median)) {
cat("\nMedian survival\n")
z <- x$median
z$estimate <- .r4vn_surv_ci_text(z$median, z$lower, z$upper, 2)
print(z[, c("group", "estimate"), drop = FALSE], row.names = FALSE)
}
if (!is.null(x$lifetable) && nrow(x$lifetable)) {
cat("\nLife table", if (isTRUE(x$metadata$competing))
" (Aalen-Johansen event history)" else " (Kaplan-Meier)", "\n", sep = "")
print(x$lifetable, row.names = FALSE)
}
if (!is.null(x$at) && nrow(x$at) && is.null(x$risk)) {
cat("\nSurvival probability at requested times\n")
z <- x$at
z$estimate <- .r4vn_surv_ci_text(z$surv, z$lower, z$upper, 2, percent = TRUE)
print(z[, intersect(c("group", "time", "n.risk", "estimate"), names(z)), drop = FALSE], row.names = FALSE)
}
if (!is.null(x$risk) && nrow(x$risk)) {
cat("\nCumulative risk", if (isTRUE(x$metadata$competing)) " (Aalen-Johansen CIF)" else " (1-KM)", "\n", sep = "")
z <- x$risk
z$estimate <- .r4vn_surv_ci_text(z$risk, z$risk_lower, z$risk_upper, 2, percent = TRUE)
print(z[, intersect(c("group", "time", "n.risk", "estimate"), names(z)), drop = FALSE], row.names = FALSE)
}
if (!is.null(x$rate) && nrow(x$rate)) {
cat("\nIncidence rate\n")
z <- x$rate
z$estimate <- .r4vn_surv_ci_text(z$rate, z$lower, z$upper, 2)
print(z[, c("group", "interval", "events", "person_time", "estimate"), drop = FALSE], row.names = FALSE)
}
if (!is.null(x$logrank)) cat("\nLog-rank: chi-square = ", .r4vn_surv_fmt(x$logrank$chisq[1], 2),
", df = ", x$logrank$df[1], ", p ", .r4vn_surv_fmt_p(x$logrank$p[1], 3), "\n", sep = "")
if (!is.null(x$risk_compare)) {
cat("\nRisk comparison\n"); print(x$risk_compare, row.names = FALSE)
}
if (!is.null(x$irr)) {
cat("\nIncidence rate ratio\n"); print(x$irr, row.names = FALSE)
}
print_model <- function(z, heading, effect = "HR") {
if (is.null(z) || !nrow(z)) return(invisible(NULL))
cat("\n", heading, "\n", sep = "")
zz <- z
zz[[paste0(effect, " (95% CI)")]] <- ifelse(zz$reference, "Ref.", .r4vn_surv_ci_text(zz$estimate, zz$lower, zz$upper, 2))
zz$p <- vapply(zz$p, .r4vn_surv_fmt_p, character(1), digits = 3)
print(zz[, intersect(c("variable", "variable_label", "level",
paste0(effect, " (95% CI)"), "p"), names(zz)), drop = FALSE],
row.names = FALSE)
}
print_model(x$cox$crude, "Crude Cox proportional hazards", "HR")
print_model(x$cox$adjusted, "Adjusted Cox proportional hazards", "HR")
print_model(x$cox$multi, "Multivariable Cox proportional hazards", "HR")
if (!is.null(x$cox$interaction) && nrow(x$cox$interaction)) {
cat("\nInteraction terms\n")
z <- x$cox$interaction
z$`HR (95% CI)` <- .r4vn_surv_ci_text(z$estimate, z$lower, z$upper, 2)
z$p <- vapply(z$p, .r4vn_surv_fmt_p, character(1), digits = 3)
print(z[, c("term", "HR (95% CI)", "p"), drop = FALSE], row.names = FALSE)
}
if (!is.null(x$cox$diagnostics) && nrow(x$cox$diagnostics)) {
z <- x$cox$diagnostics[1L, ]
cat(sprintf("\nModel N = %d; events = %d; concordance = %.3f; LR p %s\n",
as.integer(z$n), as.integer(z$events), z$concordance, .r4vn_surv_fmt_p(z$LR_p, 3)))
}
if (!is.null(x$finegray)) print_model(x$finegray$table, "Fine-Gray subdistribution hazards", "SHR")
if (!is.null(x$ph)) {
cat("\nProportional-hazards test\n")
z <- x$ph
if ("p" %in% names(z)) z$p <- vapply(z$p, .r4vn_surv_fmt_p, character(1), digits = 3)
print(z, row.names = FALSE)
}
if (!is.null(x$rmst)) {
cat("\nRestricted mean survival time\n")
print(x$rmst$table, row.names = FALSE)
if (!is.null(x$rmst$difference)) print(x$rmst$difference, row.names = FALSE)
}
if (!is.null(x$interpretation) && nrow(x$interpretation)) {
cat("\nInterpretation\n")
for (i in seq_len(nrow(x$interpretation))) {
cat("- ", x$interpretation$section[i], ": ",
x$interpretation$interpretation[i], "\n", sep = "")
}
}
invisible(x)
}
#' Direct Cox Proportional Hazards Model
#'
#' A compact R4VN wrapper around `tabsurv()` for a final multivariable Cox model.
#'
#' @inheritParams tabsurv
#' @param ph Logical; test the proportional-hazards assumption. Default `FALSE` in `cox()`. Setting `diagnosis = TRUE` also requests this diagnostic.
#' @param diagnosis Logical; if `TRUE`, append Cox-model diagnostics including concordance, proportional-hazards testing, and residual summaries. Default `FALSE`.
#' @return An object of class `r4vn_surv`.
#' @examples
#' if (requireNamespace("survival", quietly = TRUE)) {
#' d <- data.frame(
#' time = c(5, 8, 10, 12, 15, 18, 20, 22, 25, 30),
#' event = c(1, 0, 1, 1, 0, 1, 0, 1, 1, 0),
#' age = c(40, 45, 50, 55, 60, 48, 52, 63, 58, 67),
#' sex = factor(rep(c("Female", "Male"), 5))
#' )
#' cox(time, event, vars = vars(c.age, i.sex), data = d, show = FALSE)
#' cox(time, event, vars = vars(c.age, i.sex), data = d,
#' diagnosis = TRUE, show = FALSE)
#' cox(time, event, vars = vars(c.age), data = d, ph = TRUE, show = FALSE)
#' }
#' @export
cox <- function(time, event, vars, data = NULL, failure = NULL, id = NULL,
start = NULL, strata = NULL, cluster = NULL, frailty = NULL,
ph = FALSE, ci = .95, ties = c("efron", "breslow", "exact"),
diagnosis = FALSE, show = TRUE, console = FALSE) {
call <- match.call()
# Resolve the default choice before forwarding to tabsurv().
ties <- match.arg(ties)
d <- .r4vn_surv_data(data)
tn <- .r4vn_surv_name(substitute(time), d, "time")
en <- .r4vn_surv_name(substitute(event), d, "event")
idn <- .r4vn_surv_name(substitute(id), d, "id", TRUE)
stn <- .r4vn_surv_name(substitute(start), d, "start", TRUE)
srn <- .r4vn_surv_name(substitute(strata), d, "strata", TRUE)
cln <- .r4vn_surv_name(substitute(cluster), d, "cluster", TRUE)
frn <- .r4vn_surv_name(substitute(frailty), d, "frailty", TRUE)
# tabsurv() uses NSE for time/event/id/start/strata/cluster/frailty. A
# direct call such as `id = idn` would therefore expose the symbol `idn` to
# tabsurv(), even when its value is NULL, and tabsurv() would look for a
# column literally named "idn". Build the call from evaluated values so
# character variable names and literal NULLs arrive intact.
args <- list(
time = tn, event = en, vars = vars, data = d, failure = failure,
id = idn, start = stn, strata = srn, cluster = cln, frailty = frn,
km = FALSE, cox = TRUE, multi = TRUE, ph = isTRUE(ph) || isTRUE(diagnosis), ci = ci,
ties = ties, report = "custom", plot = FALSE,
interpretation = FALSE, show = FALSE, console = FALSE
)
out <- do.call(tabsurv, args)
if (!isTRUE(diagnosis) && is.list(out$cox)) out$cox$diagnostics <- NULL
if (isTRUE(diagnosis) && !is.null(out$cox$multi_fit)) {
out$model_diagnostics <- .r4vn_model_diagnosis(out$cox$multi_fit, kind = "cox")
} else if (isTRUE(diagnosis) && !is.null(out$raw$model)) {
out$model_diagnostics <- .r4vn_model_diagnosis(out$raw$model, kind = "cox")
}
out$call <- call
.r4vn_show(out, show = show, console = console)
}
# ============================================================================
# Viewer renderers for survival commands
# ============================================================================
.r4vn_surv_model_display <- function(z, effect = "HR", digit = 2, p_digit = 3) {
if (is.null(z) || !nrow(z)) return(NULL)
zz <- z
zz[[paste0(effect, " (95% CI)")]] <- ifelse(zz$reference, "Ref.",
.r4vn_surv_ci_text(zz$estimate, zz$lower, zz$upper, digit))
zz$p <- vapply(zz$p, .r4vn_surv_fmt_p, character(1), digits = p_digit)
zz[, intersect(c("variable", "variable_label", "level", "term",
paste0(effect, " (95% CI)"), "p"), names(zz)), drop = FALSE]
}
.r4vn_surv_viewer <- function(x, subtitle = "Survival analysis") {
blocks <- character()
ov <- x$overview
if (!is.null(ov) && nrow(ov)) {
overview <- data.frame(Statistic = names(ov), Value = as.character(ov[1, ]), stringsAsFactors = FALSE)
blocks <- c(blocks, paste0('<section class="r4vn-section"><h2>Overview</h2>',
.r4vn_view_key_values(overview), '</section>'))
}
if (!is.null(x$descriptive) && nrow(x$descriptive)) {
blocks <- c(blocks, .r4vn_view_section(
paste0("Outcome by ", x$metadata$by_label %||% "group"), x$descriptive
))
}
if (!is.null(x$followup) && nrow(x$followup)) {
z <- x$followup; z$`Median follow-up (95% CI)` <- .r4vn_surv_ci_text(z$median, z$lower, z$upper, 2)
blocks <- c(blocks, .r4vn_view_section("Follow-up", z[, intersect(c("group", "Median follow-up (95% CI)"), names(z)), drop = FALSE]))
}
if (!is.null(x$median) && nrow(x$median)) {
z <- x$median; z$`Median survival (95% CI)` <- .r4vn_surv_ci_text(z$median, z$lower, z$upper, 2)
blocks <- c(blocks, .r4vn_view_section("Median survival", z[, intersect(c("group", "Median survival (95% CI)"), names(z)), drop = FALSE]))
}
if (!is.null(x$lifetable) && nrow(x$lifetable)) {
ttl <- if (isTRUE(x$metadata$competing)) {
"Life table (Aalen-Johansen event history)"
} else "Life table (Kaplan-Meier)"
blocks <- c(blocks, .r4vn_view_section(ttl, x$lifetable))
}
if (!is.null(x$at) && nrow(x$at) && is.null(x$risk)) {
z <- x$at; z$`Survival (95% CI)` <- .r4vn_surv_ci_text(z$surv, z$lower, z$upper, 2, percent = TRUE)
blocks <- c(blocks, .r4vn_view_section("Survival probability at requested times", z[, intersect(c("group", "time", "n.risk", "Survival (95% CI)"), names(z)), drop = FALSE]))
}
if (!is.null(x$risk) && nrow(x$risk)) {
z <- x$risk; z$`Risk (95% CI)` <- .r4vn_surv_ci_text(z$risk, z$risk_lower, z$risk_upper, 2, percent = TRUE)
ttl <- if (isTRUE(x$metadata$competing)) "Cumulative incidence (Aalen-Johansen CIF)" else "Cumulative risk (1 - KM)"
blocks <- c(blocks, .r4vn_view_section(ttl, z[, intersect(c("group", "time", "n.risk", "Risk (95% CI)"), names(z)), drop = FALSE]))
}
if (!is.null(x$rate) && nrow(x$rate)) {
z <- x$rate; z$`Rate (95% CI)` <- .r4vn_surv_ci_text(z$rate, z$lower, z$upper, 2)
blocks <- c(blocks, .r4vn_view_section("Incidence rate", z[, intersect(c("group", "interval", "events", "person_time", "Rate (95% CI)"), names(z)), drop = FALSE]))
}
if (!is.null(x$logrank) && nrow(x$logrank)) {
z <- data.frame(
`Chi-square` = .r4vn_surv_fmt(x$logrank$chisq, 2),
df = .r4vn_surv_fmt(x$logrank$df, 0),
p = vapply(x$logrank$p, .r4vn_surv_fmt_p, character(1), digits = 3),
stringsAsFactors = FALSE, check.names = FALSE
)
blocks <- c(blocks, .r4vn_view_section("Log-rank test", z))
}
if (!is.null(x$risk_compare) && nrow(x$risk_compare)) {
z <- x$risk_compare
out <- z[, intersect(c("time", "reference", "comparison"), names(z)), drop = FALSE]
if ("rr" %in% names(z)) {
out$`RR (95% CI)` <- .r4vn_surv_ci_text(z$rr, z$rr_lower, z$rr_upper, 2)
out$`RR p` <- vapply(z$rr_p, .r4vn_surv_fmt_p, character(1), digits = 3)
}
if ("rd" %in% names(z)) {
out$`RD (95% CI)` <- .r4vn_surv_ci_text(z$rd, z$rd_lower, z$rd_upper, 3)
out$`RD p` <- vapply(z$rd_p, .r4vn_surv_fmt_p, character(1), digits = 3)
}
blocks <- c(blocks, .r4vn_view_section("Risk comparison", out))
}
if (!is.null(x$irr) && nrow(x$irr)) {
z <- x$irr
out <- z[, intersect(c("interval", "reference", "comparison"), names(z)), drop = FALSE]
out$`IRR (95% CI)` <- .r4vn_surv_ci_text(z$irr, z$lower, z$upper, 2)
out$p <- vapply(z$p, .r4vn_surv_fmt_p, character(1), digits = 3)
blocks <- c(blocks, .r4vn_view_section("Incidence rate ratio", out))
}
if (!is.null(x$cox$crude)) blocks <- c(blocks, .r4vn_view_section("Crude Cox proportional hazards", .r4vn_surv_model_display(x$cox$crude, "HR")))
if (!is.null(x$cox$adjusted)) blocks <- c(blocks, .r4vn_view_section("Adjusted Cox proportional hazards", .r4vn_surv_model_display(x$cox$adjusted, "HR")))
if (!is.null(x$cox$multi)) blocks <- c(blocks, .r4vn_view_section("Multivariable Cox proportional hazards", .r4vn_surv_model_display(x$cox$multi, "HR")))
if (!is.null(x$cox$interaction) && nrow(x$cox$interaction)) {
z <- x$cox$interaction; z$`HR (95% CI)` <- .r4vn_surv_ci_text(z$estimate, z$lower, z$upper, 2); z$p <- vapply(z$p, .r4vn_surv_fmt_p, character(1), digits = 3)
blocks <- c(blocks, .r4vn_view_section("Interaction terms", z[, intersect(c("term", "HR (95% CI)", "p"), names(z)), drop = FALSE]))
}
if (!is.null(x$cox$diagnostics) && nrow(x$cox$diagnostics)) {
z <- x$cox$diagnostics
pcols <- grep("_p$|^p$", names(z), ignore.case = TRUE, value = TRUE)
for (nm in pcols) z[[nm]] <- vapply(z[[nm]], .r4vn_surv_fmt_p, character(1), digits = 3)
for (nm in setdiff(names(z), pcols)) {
if (is.numeric(z[[nm]])) z[[nm]] <- if (nm %in% c("n", "events", "LR_df", "Wald_df", "Score_df")) .r4vn_surv_fmt(z[[nm]], 0) else .r4vn_surv_fmt(z[[nm]], 3)
}
blocks <- c(blocks, .r4vn_view_section("Cox model diagnostics", z))
}
if (!is.null(x$finegray$table)) blocks <- c(blocks, .r4vn_view_section("Fine-Gray subdistribution hazards", .r4vn_surv_model_display(x$finegray$table, "SHR")))
if (!is.null(x$model_diagnostics) && length(x$model_diagnostics)) {
for (nm in names(x$model_diagnostics)) {
z <- x$model_diagnostics[[nm]]
if (is.data.frame(z) && nrow(z)) blocks <- c(blocks, .r4vn_view_section(nm, z))
}
}
if (!is.null(x$ph) && nrow(x$ph)) {
z <- x$ph
for (nm in names(z)) {
if (tolower(nm) == "p") z[[nm]] <- vapply(z[[nm]], .r4vn_surv_fmt_p, character(1), digits = 3)
else if (is.numeric(z[[nm]])) z[[nm]] <- .r4vn_surv_fmt(z[[nm]], 3)
}
blocks <- c(blocks, .r4vn_view_section("Proportional-hazards test", z))
}
if (!is.null(x$rmst$table) && nrow(x$rmst$table)) {
z <- x$rmst$table
out <- z[, intersect(c("group", "tau"), names(z)), drop = FALSE]
if ("tau" %in% names(out)) out$tau <- .r4vn_surv_fmt(out$tau, 2)
out$`RMST (95% CI)` <- .r4vn_surv_ci_text(z$rmst, z$lower, z$upper, 2)
out$SE <- .r4vn_surv_fmt(z$se, 2)
blocks <- c(blocks, .r4vn_view_section("Restricted mean survival time", out))
}
if (!is.null(x$rmst$difference) && nrow(x$rmst$difference)) {
z <- x$rmst$difference
out <- z[, intersect(c("reference", "comparison"), names(z)), drop = FALSE]
out$`Difference (95% CI)` <- .r4vn_surv_ci_text(z$difference, z$lower, z$upper, 2)
out$p <- vapply(z$p, .r4vn_surv_fmt_p, character(1), digits = 3)
blocks <- c(blocks, .r4vn_view_section("RMST difference", out))
}
if (!is.null(x$interpretation) && nrow(x$interpretation)) {
blocks <- c(blocks, .r4vn_view_section("Interpretation", x$interpretation))
}
notes <- character()
if (isTRUE(x$metadata$missing) || x$metadata$excluded_survival > 0L)
notes <- c(notes, paste0("Excluded from survival summaries because required fields were incomplete: ", x$metadata$excluded_survival, "."))
if (!is.null(x$metadata$excluded_model_base) && x$metadata$excluded_model_base != x$metadata$excluded_survival)
notes <- c(notes, paste0("Excluded from the structural model base: ", x$metadata$excluded_model_base, "."))
.r4vn_view_document(x$title %||% "R4VN survival analysis", paste0(blocks, collapse = ""), notes = notes,
subtitle = subtitle, prefix = "r4vn-tabsurv-")
}
.r4vn_surv_hierarchical_viewer <- function(x) {
results <- x$hierarchical_results
labels <- names(results) %||% x$strata_labels %||% paste0("Stratum ", seq_along(results))
summary_rows <- lapply(seq_along(results), function(i) {
ov <- results[[i]]$overview
data.frame(
Stratum = labels[i], N = ov$N[1L], Events = ov$Events[1L],
Censored = ov$Censored[1L], `Event %` = .r4vn_surv_fmt(ov$`Event percent`[1L], 1),
`Person-time` = .r4vn_surv_fmt(ov$`Person-time`[1L], 2),
stringsAsFactors = FALSE, check.names = FALSE
)
})
blocks <- .r4vn_view_section("Hierarchy overview", do.call(rbind, summary_rows))
for (i in seq_along(results)) {
z <- results[[i]]
blocks <- paste0(
blocks, '<section class="r4vn-section"><h2>',
.r4vn_view_escape(labels[i]), '</h2></section>'
)
if (!is.null(z$descriptive) && nrow(z$descriptive)) {
blocks <- paste0(blocks, .r4vn_view_section("Outcome by group", z$descriptive))
}
if (!is.null(z$median) && nrow(z$median)) {
q <- z$median
q$`Median survival (95% CI)` <- .r4vn_surv_ci_text(q$median, q$lower, q$upper, 2)
blocks <- paste0(blocks, .r4vn_view_section(
"Median survival", q[, intersect(c("group", "Median survival (95% CI)"), names(q)), drop = FALSE]
))
}
if (!is.null(z$lifetable) && nrow(z$lifetable)) {
blocks <- paste0(blocks, .r4vn_view_section(
if (isTRUE(z$metadata$competing)) {
"Life table (Aalen-Johansen event history)"
} else "Life table (Kaplan-Meier)",
z$lifetable
))
}
if (!is.null(z$risk) && nrow(z$risk)) {
q <- z$risk
q$`Risk (95% CI)` <- .r4vn_surv_ci_text(q$risk, q$risk_lower, q$risk_upper, 2, percent = TRUE)
blocks <- paste0(blocks, .r4vn_view_section(
if (isTRUE(z$metadata$competing)) "Cumulative incidence" else "Cumulative risk",
q[, intersect(c("group", "time", "n.risk", "Risk (95% CI)"), names(q)), drop = FALSE]
))
}
if (!is.null(z$rate) && nrow(z$rate)) {
q <- z$rate
q$`Rate (95% CI)` <- .r4vn_surv_ci_text(q$rate, q$lower, q$upper, 2)
blocks <- paste0(blocks, .r4vn_view_section(
"Incidence rate", q[, intersect(c("group", "interval", "events", "person_time", "Rate (95% CI)"), names(q)), drop = FALSE]
))
}
if (!is.null(z$logrank) && nrow(z$logrank)) {
q <- data.frame(
`Chi-square` = .r4vn_surv_fmt(z$logrank$chisq, 2),
df = .r4vn_surv_fmt(z$logrank$df, 0),
p = vapply(z$logrank$p, .r4vn_surv_fmt_p, character(1), digits = 3),
stringsAsFactors = FALSE, check.names = FALSE
)
blocks <- paste0(blocks, .r4vn_view_section("Log-rank test", q))
}
model <- if (!is.null(z$cox$multi)) z$cox$multi else z$cox$crude
if (!is.null(model) && nrow(model)) {
blocks <- paste0(blocks, .r4vn_view_section(
if (!is.null(z$cox$multi)) "Multivariable Cox proportional hazards" else "Crude Cox proportional hazards",
.r4vn_surv_model_display(model, "HR")
))
}
if (!is.null(z$interpretation) && nrow(z$interpretation)) {
blocks <- paste0(blocks, .r4vn_view_section("Interpretation", z$interpretation))
}
}
notes <- paste0(
"Hierarchical by: ", paste(x$hierarchical_by$strata_labels %||%
x$hierarchical_by$strata, collapse = " > "),
"; innermost survival group: ", x$hierarchical_by$by_label %||%
x$hierarchical_by$by, "."
)
.r4vn_view_document(
x$title %||% "R4VN hierarchical survival analysis", blocks,
notes = notes, subtitle = "Hierarchical survival analysis",
prefix = "r4vn-tabsurv-hierarchical-"
)
}
.r4vn_viewer_tabsurv <- function(x) {
if (!is.null(x$hierarchical_results)) return(.r4vn_surv_hierarchical_viewer(x))
.r4vn_surv_viewer(x, "Survival, incidence and time-to-event analysis")
}
.r4vn_viewer_cox <- function(x) .r4vn_surv_viewer(x, "Cox proportional hazards model")
Try the R4VN package in your browser
Any scripts or data that you put into this service are public.
R4VN documentation built on Sept. 30, 2026, 5:13 p.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Defines functions print.r4vn_surv .r4vn_tabsurv_plot .r4vn_tabsurv_contract .r4vn_tabsurv_hierarchical_result .r4vn_tabsurv_profile .r4vn_tabsurv_by_info .r4vn_tabsurv_core tabsurv
Documented in tabsurv
#' Comprehensive Survival Analysis Table
#'
#' Performs descriptive survival analysis, Kaplan-Meier/Aalen-Johansen estimates,
#' optional life tables, cumulative incidence at selected times, incidence rate,
#' log-rank tests, Cox regression, proportional-hazards diagnostics, RMST,
#' competing-risk Fine-Gray models, and counting-process/recurrent-event Cox models.
#'
#' @param time Follow-up or stop-time variable, supplied without quotes.
#' @param event Event/status variable, supplied without quotes.
#' @param vars Optional predictor specification created by `vars()`.
#' @param by Optional grouping variable for survival curves and comparisons.
#' Hierarchical syntax is supported: in `by = vars(province, sex, treatment)`,
#' treatment is the innermost curve/comparison group and province > sex are
#' ordered outer strata.
#' @param data Optional data frame. When omitted, active R4VN data are used.
#' @param failure Value of `event` representing the event of interest. For a
#' binary event it defaults to the second factor level or larger numeric value.
#' @param compete Optional competing-event value(s). When supplied, `risk = TRUE`
#' uses the Aalen-Johansen cumulative incidence function.
#' @param id Optional subject identifier for counting-process/recurrent data.
#' @param start Optional start/entry time. When supplied, `time` is treated as stop time.
#' @param unit Optional display unit such as "day", "month", or "year".
#' @param followup Estimate median follow-up using reverse Kaplan-Meier when possible.
#' With `report = "auto"`, this is enabled unless explicitly set to `FALSE`.
#' @param km Fit Kaplan-Meier (ordinary survival) or Aalen-Johansen (competing risks).
#' @param lifetable Show a detailed life table at every observed time. The
#' default is `FALSE`. For ordinary survival, the table reports numbers at
#' risk, events, censoring, conditional survival, cumulative Kaplan-Meier
#' survival, cumulative risk, standard error, and confidence limits. With
#' competing risks, it reports the corresponding Aalen-Johansen event-history
#' table and cumulative incidence.
#' @param at Optional time points for survival/risk/rate summaries. With
#' `report = "auto"` or `"full"`, three representative follow-up times are
#' selected automatically when `at` is omitted.
#' @param risk Report cumulative risk at `at`. For ordinary survival this is 1-S(t);
#' with competing risks it is the cumulative incidence function.
#' @param cuminc Optional numeric time points at which cumulative incidence is
#' required, for example `cuminc = c(6, 12, 24)`. This directly activates
#' cumulative-risk output without also requiring `risk = TRUE`. Ordinary
#' survival uses 1-KM; competing-risk analysis uses the Aalen-Johansen CIF.
#' @param rate `FALSE`, `TRUE`, `"overall"`, `"cumulative"`, `"interval"`,
#' or `"all"`. `"all"` reports overall, cumulative, and interval-specific
#' rates. The automatic profile uses the overall incidence rate.
#' @param scale Rate multiplier, e.g. 100 for events per 100 person-time units.
#' @param logrank Perform a log-rank test when `by` is supplied and no competing risk exists.
#' @param rr,rd Compare cumulative risks between two `by` groups using approximate
#' risk ratio or risk difference inference based on survival-estimate standard errors.
#' @param irr Compare incidence rates between two `by` groups.
#' @param cox Fit crude Cox models for variables in `vars`.
#' @param adjusted FALSE/NULL, TRUE (adjust each focal predictor for all other focal
#' predictors), or a `vars()`/character set of adjustment covariates.
#' @param multi FALSE/NULL, TRUE (all `vars` in one model), or a `vars()`/character
#' set defining the final multivariable Cox model.
#' @param strata Optional stratification variable for Cox regression.
#' @param cluster Optional clustering variable for robust Cox variance.
#' @param frailty Optional shared-frailty variable. Do not combine with `cluster`.
#' @param finegray Fit a Fine-Gray subdistribution hazards model when `compete` is supplied.
#' @param recurrent FALSE/TRUE or "ag". TRUE is Andersen-Gill and requires
#' `id` and `start`. Automatic profiles suppress ordinary KM/RMST modules for
#' recurrent-event data unless the user explicitly requests them.
#' @param rmst Compute restricted mean survival time.
#' @param tau Restriction time for RMST. Defaults to the largest common curve time.
#' @param ph Test the proportional-hazards assumption with `cox.zph()` for the final Cox model.
#' @param interaction Optional `vars(a, b)` containing exactly two variables to include
#' their interaction in the final multivariable Cox model.
#' @param superby Optional outer subgroup variable retained for backward compatibility.
#' For new code, multiple ordered outer strata can be supplied directly in
#' `by = vars(stratum1, stratum2, group)`.
#' @param ci Confidence level, default 0.95.
#' @param digit,p_digit,effect_digit Display digits.
#' @param missing Show missing/exclusion information when printing.
#' @param plot Draw a survival/CIF curve using `gsurv()` after analysis. In the
#' automatic profiles, the plot includes confidence limits, the log-rank
#' p-value when available, and a number-at-risk table.
#' @param title Optional title.
#' @param show Logical; open the formatted result in the Viewer. Default `TRUE`.
#' @param console Logical; also print the traditional result in the Console. Default `FALSE`.
#' @param ai Prepare a compact de-identified interpretation payload in `$ai_text`.
#' @param ties Cox tie method: "efron", "breslow", or "exact".
#' @param report Reporting profile: `"auto"` (context-sensitive comprehensive
#' output), `"brief"` (descriptive survival summary), `"full"` (all valid
#' modules), or `"custom"` (backward-compatible concise defaults plus explicitly requested
#' modules).
#' @param plot_args Named list of additional arguments passed to `gsurv()`.
#' @param interpretation Add a cautious, deterministic interpretation table.
#' The default is `FALSE`; use `interpretation = TRUE` when narrative output
#' is wanted.
#' @param export Optional export format accepted by `tabexport()`, such as
#' `"docx"`, `"xlsx"`, or `"html"`.
#' @param file Optional export filename. Its extension may also determine the
#' export format.
#' @param open Open the exported file when supported.
#' @param strict If `TRUE`, an unavailable optional module stops the analysis.
#' The default `FALSE` keeps the main report and records a warning instead.
#'
#' @return An object of class `r4vn_surv`. Backward-compatible components are
#' retained, with a consistent reporting contract in `$descriptive`,
#' `$estimates`, `$tests`, `$diagnostics`, `$interpretation`, `$tables`,
#' `$plots`, `$models`, `$metadata`, and `$call`.
#' @export
#'
#' @examples
#' if (requireNamespace("survival", quietly = TRUE)) {
#' # Reproducible two-group data from the survival package.
#' d <- survival::lung
#' d$death <- as.integer(d$status == 2)
#' d$group <- factor(d$sex, levels = c(1, 2),
#' labels = c("Male", "Female"))
#'
#' # 1. Complete two-group report. This includes the log-rank test.
#' km <- tabsurv(
#' time, death, by = group, data = d, failure = 1,
#' unit = "day", at = c(90, 180, 365, 540),
#' report = "auto", plot = FALSE, show = FALSE
#' )
#' km$logrank
#' km$tests$logrank
#' km$logrank$p
#'
#' \donttest{
#' # 2. Cumulative incidence at 6, 12, and 24 months.
#' d$month <- d$time / 30.4375
#' ci_month <- tabsurv(
#' month, death, by = group, data = d, failure = 1,
#' cuminc = c(6, 12, 24), report = "custom", show = FALSE
#' )
#' ci_month$cuminc
#'
#' # 3. Detailed life table and interpretation are both opt-in.
#' km_detail <- tabsurv(
#' time, death, by = group, data = d, failure = 1,
#' at = c(90, 180, 365, 540),
#' lifetable = TRUE, interpretation = TRUE, show = FALSE
#' )
#' head(km_detail$lifetable)
#' km_detail$interpretation
#'
#' # 4. A compact KM plus log-rank analysis without automatic extras.
#' km_simple <- tabsurv(
#' time, death, by = group, data = d, failure = 1,
#' report = "custom", km = TRUE, logrank = TRUE,
#' plot = FALSE, show = FALSE
#' )
#'
#' # 5. Explicit two-group effect measures and RMST.
#' km_compare <- tabsurv(
#' time, death, by = group, data = d, failure = 1,
#' at = c(90, 180, 365, 540), risk = TRUE,
#' rr = TRUE, rd = TRUE, rate = "all", irr = TRUE,
#' rmst = TRUE, tau = 365, show = FALSE
#' )
#' km_compare$risk_compare
#' km_compare$irr
#' km_compare$rmst
#'
#' # 6. Publication graphs, including risk tables, are documented in ?gsurv.
#' # Keeping graphics out of this example also keeps tabsurv() examples fast
#' # and executable on non-interactive CRAN check devices.
#'
#' # 7. With competing risks, use the Aalen-Johansen CIF, not 1-KM.
#' set.seed(2026)
#' n <- 180
#' t1 <- rexp(n, 0.07)
#' t2 <- rexp(n, 0.05)
#' tc <- runif(n, 4, 30)
#' tm <- pmin(t1, t2, tc)
#' dcr <- data.frame(
#' time = tm,
#' status = ifelse(tm == t1, 1L, ifelse(tm == t2, 2L, 0L)),
#' group = factor(rep(c("A", "B"), each = n / 2))
#' )
#' cif <- tabsurv(
#' time, status, by = group, data = dcr,
#' failure = 1, compete = 2, cuminc = c(6, 12, 24),
#' report = "custom", show = FALSE
#' )
#' cif$cuminc
#' # See ?gsurv for publication CIF graphs and risk tables.
#' }
#' }
tabsurv <- function(time, event, vars = NULL, by = NULL, data = NULL,
failure = NULL, compete = NULL,
id = NULL, start = NULL, unit = NULL,
followup = NULL, km = NULL, lifetable = FALSE, at = NULL,
risk = NULL, cuminc = NULL, rate = NULL, scale = 100,
logrank = NULL, rr = NULL, rd = NULL, irr = NULL,
cox = NULL, adjusted = NULL, multi = NULL,
strata = NULL, cluster = NULL, frailty = NULL,
finegray = NULL, recurrent = FALSE,
rmst = NULL, tau = NULL,
ph = NULL, interaction = FALSE, superby = NULL,
ci = 0.95, digit = 2, p_digit = 3, effect_digit = 2,
missing = FALSE, plot = NULL,
title = NULL, show = TRUE, console = FALSE, ai = FALSE,
ties = c("efron", "breslow", "exact"),
report = c("auto", "brief", "full", "custom"),
plot_args = list(), interpretation = FALSE,
export = NULL, file = NULL, open = FALSE,
strict = FALSE) {
.r4vn_surv_require()
call <- match.call()
env <- parent.frame()
d <- .r4vn_surv_data(data)
report <- match.arg(report)
ties <- match.arg(ties)
if (!is.list(plot_args) ||
(length(plot_args) && (is.null(names(plot_args)) || any(!nzchar(names(plot_args)))))) {
stop("`plot_args` must be a named list.", call. = FALSE)
}
if (!is.logical(interpretation) || length(interpretation) != 1L || is.na(interpretation)) {
stop("`interpretation` must be TRUE or FALSE.", call. = FALSE)
}
if (!is.logical(lifetable) || length(lifetable) != 1L || is.na(lifetable)) {
stop("`lifetable` must be TRUE or FALSE.", call. = FALSE)
}
if (!is.null(cuminc) && (!is.numeric(cuminc) || !length(cuminc) ||
any(!is.finite(cuminc)) || any(cuminc < 0))) {
stop("`cuminc` must contain non-negative finite time points.", call. = FALSE)
}
if (!is.null(cuminc)) cuminc <- sort(unique(cuminc))
if (!is.logical(strict) || length(strict) != 1L || is.na(strict)) {
stop("`strict` must be TRUE or FALSE.", call. = FALSE)
}
by_expr <- substitute(by)
by_info <- .r4vn_tabsurv_by_info(by_expr, d, env)
time_name <- .r4vn_surv_name(substitute(time), d, "time")
provided <- c(
followup = !missing(followup), km = !missing(km), at = !missing(at),
risk = !missing(risk), rate = !missing(rate), logrank = !missing(logrank),
rr = !missing(rr), rd = !missing(rd), irr = !missing(irr),
cox = !missing(cox), adjusted = !missing(adjusted), multi = !missing(multi),
finegray = !missing(finegray), rmst = !missing(rmst), ph = !missing(ph),
plot = !missing(plot)
)
profile <- .r4vn_tabsurv_profile(
report = report, provided = provided,
values = list(
followup = followup, km = km, at = at, risk = risk, rate = rate,
logrank = logrank, rr = rr, rd = rd, irr = irr, cox = cox,
adjusted = adjusted, multi = multi, finegray = finegray, rmst = rmst,
ph = ph, plot = plot, vars = vars
),
data = d, time_name = time_name, by_name = by_info$by,
competing = !is.null(compete), recurrent = recurrent
)
for (nm in names(profile)) assign(nm, profile[[nm]])
if (!is.null(cuminc)) risk <- TRUE
core_call <- call
core_call[[1L]] <- quote(.r4vn_tabsurv_core)
core_call$data <- quote(.r4vn_tabsurv_data)
core_call$followup <- followup
core_call$km <- km
core_call$lifetable <- lifetable
core_call$at <- at
core_call$risk <- risk
core_call$cuminc <- cuminc
core_call$rate <- rate
core_call$logrank <- logrank
core_call$rr <- rr
core_call$rd <- rd
core_call$irr <- irr
core_call$cox <- cox
if (!isTRUE(provided[["adjusted"]]) || is.logical(adjusted)) core_call$adjusted <- adjusted
if (!isTRUE(provided[["multi"]]) || is.logical(multi)) core_call$multi <- multi
core_call$finegray <- finegray
core_call$rmst <- rmst
core_call$ph <- ph
core_call$plot <- plot
core_call$show <- FALSE
core_call$console <- FALSE
core_call$ties <- ties
core_call$report <- report
core_call$strict <- strict
core_call$plot_args <- NULL
core_call$interpretation <- NULL
core_call$export <- NULL
core_call$file <- NULL
core_call$open <- NULL
if (!is.null(by_info$spec)) {
spec <- by_info$spec
outer <- spec$strata
super_expr <- substitute(superby)
if (!.r4vn_expr_is_null(super_expr)) {
explicit <- .r4vn_resolve_name_spec(
super_expr, d, env, "superby", allow_null = TRUE, multiple = TRUE
)
outer <- unique(c(explicit, outer))
}
outer <- setdiff(outer, spec$by)
if (length(outer)) {
ids <- .r4vn_strata_indices(d, outer)
if (!length(ids)) stop("No complete strata are available for survival analysis.", call. = FALSE)
labels <- vapply(ids, function(idx) .r4vn_stratum_label(d, outer, idx), character(1))
core_call$by <- as.name(spec$by)
core_call$superby <- NULL
core_call$plot <- FALSE
core_call$ai <- FALSE
results <- lapply(ids, function(idx) {
ee <- new.env(parent = env)
ee$.r4vn_tabsurv_core <- .r4vn_tabsurv_core
ee$.r4vn_tabsurv_data <- d[idx, , drop = FALSE]
z <- eval(core_call, envir = ee)
z$call <- call
.r4vn_tabsurv_contract(z, interpretation = interpretation)
})
names(results) <- labels
out <- .r4vn_tabsurv_hierarchical_result(
results, labels, spec, outer, d, call, title, missing, ci, scale, report
)
if (isTRUE(ai)) {
out$ai_text <- lapply(results, function(z) {
tryCatch(.r4vn_surv_ai_text(z, effect_digit, p_digit),
error = function(e) NULL)
})
}
} else {
core_call$by <- as.name(spec$by)
core_call$superby <- NULL
ee <- new.env(parent = env)
ee$.r4vn_tabsurv_core <- .r4vn_tabsurv_core
ee$.r4vn_tabsurv_data <- d
out <- eval(core_call, envir = ee)
}
} else {
ee <- new.env(parent = env)
ee$.r4vn_tabsurv_core <- .r4vn_tabsurv_core
ee$.r4vn_tabsurv_data <- d
out <- eval(core_call, envir = ee)
}
out$call <- call
out <- .r4vn_tabsurv_contract(out, interpretation = interpretation)
out <- .r4vn_tabsurv_plot(
out, plot = plot, plot_args = plot_args, report = report, show = show
)
if (!is.null(export) || !is.null(file)) {
out$export <- survexport(
out, export = export, file = file, open = open,
title = title %||% "R4VN survival analysis"
)
}
.r4vn_show(out, show = show, console = console)
}
.r4vn_tabsurv_core <- function(time, event, vars = NULL, by = NULL, data = NULL,
failure = NULL, compete = NULL,
id = NULL, start = NULL, unit = NULL,
followup = TRUE, km = TRUE, lifetable = FALSE, at = NULL,
risk = FALSE, cuminc = NULL, rate = FALSE, scale = 100,
logrank = TRUE, rr = FALSE, rd = FALSE, irr = FALSE,
cox = FALSE, adjusted = FALSE, multi = FALSE,
strata = NULL, cluster = NULL, frailty = NULL,
finegray = FALSE, recurrent = FALSE,
rmst = FALSE, tau = NULL,
ph = FALSE, interaction = FALSE, superby = NULL,
ci = 0.95, digit = 2, p_digit = 3, effect_digit = 2,
missing = FALSE, plot = FALSE,
title = NULL, show = FALSE, console = FALSE, ai = FALSE,
ties = c("efron", "breslow", "exact"),
report = "custom", strict = FALSE) {
.r4vn_surv_require()
call <- match.call()
env <- parent.frame()
data <- .r4vn_surv_data(data)
ties <- match.arg(ties)
if (!is.numeric(ci) || length(ci) != 1L || is.na(ci) || ci <= 0 || ci >= 1) stop("`ci` must be between 0 and 1.", call. = FALSE)
if (!is.numeric(scale) || length(scale) != 1L || is.na(scale) || scale <= 0) stop("`scale` must be a positive number.", call. = FALSE)
if (!is.null(at)) {
if (!is.numeric(at) || any(!is.finite(at)) || any(at < 0)) stop("`at` must contain non-negative finite time points.", call. = FALSE)
at <- sort(unique(at))
}
if (!is.null(cuminc)) {
if (!is.numeric(cuminc) || !length(cuminc) || any(!is.finite(cuminc)) ||
any(cuminc < 0)) {
stop("`cuminc` must contain non-negative finite time points.", call. = FALSE)
}
cuminc <- sort(unique(cuminc))
risk <- TRUE
}
time_name <- .r4vn_surv_name(substitute(time), data, "time")
event_name <- .r4vn_surv_name(substitute(event), data, "event")
by_name <- .r4vn_surv_name(substitute(by), data, "by", TRUE)
id_name <- .r4vn_surv_name(substitute(id), data, "id", TRUE)
start_name <- .r4vn_surv_name(substitute(start), data, "start", TRUE)
strata_name0 <- .r4vn_surv_name(substitute(strata), data, "strata", TRUE)
cluster_name0 <- .r4vn_surv_name(substitute(cluster), data, "cluster", TRUE)
frailty_name0 <- .r4vn_surv_name(substitute(frailty), data, "frailty", TRUE)
superby_name <- .r4vn_surv_name(substitute(superby), data, "superby", TRUE)
if (!is.null(cluster_name0) && !is.null(frailty_name0)) stop("Use either `cluster` or `frailty`, not both in the same Cox model.", call. = FALSE)
if (!isFALSE(recurrent)) {
rec <- if (isTRUE(recurrent)) "ag" else tolower(as.character(recurrent)[1L])
if (!identical(rec, "ag")) stop("This implementation currently supports recurrent = TRUE or recurrent = 'ag' (Andersen-Gill).", call. = FALSE)
if (is.null(id_name) || is.null(start_name)) stop("Andersen-Gill recurrent-event analysis requires both `id` and `start`.", call. = FALSE)
if (is.null(cluster_name0)) cluster_name0 <- id_name
cox <- TRUE
multi <- if (isFALSE(multi)) TRUE else multi
}
status <- .r4vn_surv_event_status(data[[event_name]], failure, compete)
.r4vn_surv_validate_time(data[[time_name]], status$event, if (is.null(start_name)) NULL else data[[start_name]])
focal_spec <- .r4vn_surv_as_spec(vars, data, "vars", TRUE)
adjusted_spec <- NULL
if (!isFALSE(adjusted) && !is.null(adjusted)) {
adjusted_spec <- if (isTRUE(adjusted)) TRUE else .r4vn_surv_as_spec(adjusted, data, "adjusted", FALSE)
}
multi_spec <- NULL
if (!isFALSE(multi) && !is.null(multi)) {
multi_spec <- if (isTRUE(multi)) TRUE else .r4vn_surv_as_spec(multi, data, "multi", FALSE)
}
interaction_vars <- NULL
if (!isFALSE(interaction) && !is.null(interaction)) {
isp <- .r4vn_surv_as_spec(interaction, data, "interaction", FALSE)
if (nrow(isp) != 2L) stop("`interaction` must specify exactly two variables, for example `vars(treatment, sex)`.", call. = FALSE)
interaction_vars <- isp$variable
multi <- if (isFALSE(multi)) TRUE else multi
multi_spec <- if (is.null(multi_spec)) TRUE else multi_spec
}
all_spec <- .r4vn_surv_merge_specs(
focal_spec,
if (isTRUE(adjusted_spec)) NULL else adjusted_spec,
if (isTRUE(multi_spec)) NULL else multi_spec
)
if (isTRUE(adjusted_spec) || isTRUE(multi_spec)) all_spec <- .r4vn_surv_merge_specs(all_spec, focal_spec)
internal <- .r4vn_surv_internal_data(data, all_spec)
d <- internal$data
map <- internal$map
d$.time <- data[[time_name]]
d$.event <- status$event
d$.status_ms <- status$status_ms
if (!is.null(start_name)) d$.start <- data[[start_name]]
if (!is.null(by_name)) d$.by <- if (is.factor(data[[by_name]])) droplevels(data[[by_name]]) else factor(data[[by_name]], levels = unique(data[[by_name]][!is.na(data[[by_name]])]))
if (!is.null(id_name)) d$.id <- data[[id_name]]
if (!is.null(strata_name0)) d$.strata <- data[[strata_name0]]
if (!is.null(cluster_name0)) d$.cluster <- data[[cluster_name0]]
if (!is.null(frailty_name0)) d$.frailty <- data[[frailty_name0]]
if (!is.null(superby_name)) d$.superby <- if (is.factor(data[[superby_name]])) droplevels(data[[superby_name]]) else factor(data[[superby_name]])
# Use separate analysis sets. Descriptive survival estimates must not lose
# participants merely because a Cox covariate is missing.
competing <- !is.null(compete)
if (competing && !is.null(start_name) && is.null(id_name)) {
stop("Competing-risk start-stop data require `id`.", call. = FALSE)
}
desc_required <- c(".time", ".event")
if (!is.null(start_name)) desc_required <- c(desc_required, ".start")
if (!is.null(by_name) && (km || isTRUE(lifetable) || risk || !is.null(cuminc) ||
!identical(rate, FALSE) || logrank || rr || rd || irr || rmst ||
isTRUE(plot))) {
desc_required <- c(desc_required, ".by")
}
if (!is.null(start_name) && !is.null(id_name)) desc_required <- c(desc_required, ".id")
cc_desc <- .r4vn_surv_complete(d, desc_required)
d_desc <- cc_desc$data
if (!nrow(d_desc)) stop("No complete observations remain for survival-time analysis.", call. = FALSE)
if (".by" %in% names(d_desc)) d_desc$.by <- droplevels(d_desc$.by)
# Cox/Fine-Gray models use their own base set. Predictor-specific missingness
# is left to coxph()/finegray(), so crude models are not forced to use the
# complete-case sample of every other predictor.
model_required <- c(".time", ".event")
if (!is.null(start_name)) model_required <- c(model_required, ".start")
if (!is.null(strata_name0)) model_required <- c(model_required, ".strata")
if (!is.null(cluster_name0)) model_required <- c(model_required, ".cluster")
if (!is.null(frailty_name0)) model_required <- c(model_required, ".frailty")
if (!is.null(superby_name)) model_required <- c(model_required, ".superby")
if ((!isFALSE(recurrent) || (isTRUE(finegray) && !is.null(start_name))) && !is.null(id_name)) model_required <- c(model_required, ".id")
cc_model <- .r4vn_surv_complete(d, model_required)
d_model <- cc_model$data
if (!nrow(d_model)) stop("No complete observations remain for the requested model structure.", call. = FALSE)
if (!is.null(map) && nrow(map)) {
for (i in seq_len(nrow(map))) if (map$type[i] == "categorical") d_model[[map$internal[i]]] <- droplevels(d_model[[map$internal[i]]])
}
if (".superby" %in% names(d_model)) d_model$.superby <- droplevels(d_model$.superby)
overview <- .r4vn_surv_overview(d_desc)
descriptive <- .r4vn_surv_group_overview(d_desc)
followup_table <- if (isTRUE(followup)) .r4vn_surv_reverse_followup(d_desc, ci) else NULL
fit <- curve <- median_table <- at_table <- risk_table <- life_table <- NULL
logrank_table <- NULL
risk_times <- if (!is.null(cuminc)) cuminc else at
if (km || risk || rmst || isTRUE(lifetable) || !is.null(at) ||
!is.null(cuminc) || isTRUE(plot)) {
rhs <- if (".by" %in% names(d_desc)) ".by" else "1"
if (!competing) {
sf <- .r4vn_surv_formula(".start" %in% names(d_desc), rhs)
fit <- survival::survfit(sf, data = d_desc, conf.int = ci,
id = if (".start" %in% names(d_desc) && ".id" %in% names(d_desc)) d_desc$.id else NULL)
curve <- .r4vn_surv_curve_km(fit)
median_table <- .r4vn_surv_median(fit)
if (!is.null(at)) {
at_table <- .r4vn_surv_km_at(fit, at)
}
if (isTRUE(risk) && !is.null(risk_times)) {
rs <- if (!is.null(at_table) && identical(risk_times, at)) {
at_table
} else .r4vn_surv_km_at(fit, risk_times)
if (nrow(rs)) {
rs$risk <- 1 - rs$surv
rs$risk_lower <- 1 - rs$upper
rs$risk_upper <- 1 - rs$lower
rs$se <- rs$std.err
risk_table <- rs[, intersect(c("group", "time", "n.risk", "risk", "se", "risk_lower", "risk_upper", "surv", "lower", "upper"), names(rs)), drop = FALSE]
}
}
} else {
sf <- .r4vn_surv_ms_formula(".start" %in% names(d_desc), rhs)
fit <- survival::survfit(sf, data = d_desc, conf.int = ci, id = if (".id" %in% names(d_desc)) d_desc$.id else NULL)
curve <- .r4vn_surv_curve_cif(fit, "failure")
if (isTRUE(risk) && !is.null(risk_times)) {
risk_table <- .r4vn_surv_cif_at(curve, d_desc, risk_times)
names(risk_table)[names(risk_table) == "cif"] <- "risk"
# Use the Aalen-Johansen standard error; keep a CI-based fallback for
# older survival objects that do not expose std.err.
zcrit <- stats::qnorm(1 - (1 - ci) / 2)
risk_table$se <- risk_table$std.err
bad_se <- !is.finite(risk_table$se)
risk_table$se[bad_se] <- (risk_table$upper[bad_se] - risk_table$lower[bad_se]) / (2 * zcrit)
risk_table$risk_lower <- risk_table$lower
risk_table$risk_upper <- risk_table$upper
}
}
if (isTRUE(lifetable) && !is.null(curve) && nrow(curve)) {
life_table <- .r4vn_surv_lifetable(curve, competing = competing)
}
}
if (isTRUE(logrank) && ".by" %in% names(d_desc) && ".start" %in% names(d_desc) && !competing) {
warning("Log-rank testing is not reported for start-stop/counting-process data; use the Cox model for comparison.", call. = FALSE)
}
if (isTRUE(logrank) && ".by" %in% names(d_desc) && nlevels(d_desc$.by) > 1L && !competing && !(".start" %in% names(d_desc))) {
lf <- .r4vn_surv_formula(FALSE, ".by")
lr <- survival::survdiff(lf, data = d_desc)
df <- max(1L, length(lr$n) - 1L)
logrank_table <- data.frame(chisq = unname(lr$chisq), df = df,
p = stats::pchisq(lr$chisq, df, lower.tail = FALSE), stringsAsFactors = FALSE)
}
rate_table <- irr_table <- NULL
if (!identical(rate, FALSE)) {
mode <- if (is.character(rate)) tolower(rate[1L]) else TRUE
allowed_rate <- c("overall", "cumulative", "interval", "all")
if (is.character(mode) && !mode %in% allowed_rate) {
stop("`rate` must be FALSE, TRUE, 'overall', 'cumulative', 'interval', or 'all'.", call. = FALSE)
}
rate_table <- .r4vn_surv_rate_table(d_desc, at, mode, scale, ci)
if (isTRUE(irr)) {
if (!".by" %in% names(d_desc)) stop("`irr = TRUE` requires `by`.", call. = FALSE)
irr_table <- tryCatch(
.r4vn_surv_compare_rates(rate_table, ci),
error = function(e) {
if (isTRUE(strict)) stop(e)
warning("Incidence-rate comparison was not available: ", conditionMessage(e), call. = FALSE)
NULL
}
)
}
}
risk_compare <- NULL
if (isTRUE(rr) || isTRUE(rd)) {
if (is.null(risk_times) || !length(risk_times)) {
stop("`rr`/`rd` require `at` or `cuminc` time points.", call. = FALSE)
}
if (!".by" %in% names(d_desc)) stop("`rr`/`rd` require `by`.", call. = FALSE)
if (is.null(risk_table)) stop("Risk estimates are unavailable.", call. = FALSE)
risk_compare <- tryCatch(
.r4vn_surv_compare_risk(risk_table, ci, rr, rd),
error = function(e) {
if (isTRUE(strict)) stop(e)
warning("Cumulative-risk comparison was not available: ", conditionMessage(e), call. = FALSE)
NULL
}
)
}
rmst_result <- NULL
if (isTRUE(rmst)) {
if (competing) stop("RMST in this function is for ordinary survival. For competing risks use CIF/Fine-Gray results.", call. = FALSE)
if (is.null(fit)) {
rhs <- if (".by" %in% names(d_desc)) ".by" else "1"
fit <- survival::survfit(.r4vn_surv_formula(".start" %in% names(d_desc), rhs), data = d_desc, conf.int = ci)
}
if (is.null(tau)) {
if (is.null(fit$strata)) tau <- max(fit$time, na.rm = TRUE) else {
ends <- cumsum(as.integer(fit$strata)); starts <- c(1L, head(ends, -1L) + 1L)
tau <- min(vapply(Map(seq.int, starts, ends), function(ii) max(fit$time[ii], na.rm = TRUE), numeric(1)))
}
}
if (!is.numeric(tau) || length(tau) != 1L || !is.finite(tau) || tau <= 0) stop("`tau` must be one positive finite time point.", call. = FALSE)
rmst_result <- tryCatch(
.r4vn_surv_rmst(fit, tau, ci),
error = function(e) {
if (isTRUE(strict)) stop(e)
warning("RMST was not available: ", conditionMessage(e), call. = FALSE)
NULL
}
)
}
# Cox models -------------------------------------------------------------
cox_result <- list(crude = NULL, adjusted = NULL, multi = NULL,
interaction = NULL, diagnostics = NULL,
crude_fits = list(), adjusted_fits = list(), multi_fit = NULL)
need_cox <- isTRUE(cox) || !isFALSE(adjusted) || !isFALSE(multi) || isTRUE(ph) || !isFALSE(recurrent)
if (need_cox) {
if (is.null(focal_spec) || !nrow(focal_spec)) stop("Cox analysis requires `vars`.", call. = FALSE)
focal_map <- map[match(focal_spec$variable, map$variable), , drop = FALSE]
adj_obj <- if (isTRUE(adjusted_spec)) TRUE else if (is.null(adjusted_spec)) NULL else adjusted_spec
mul_obj <- if (isTRUE(multi_spec) || (isTRUE(multi) && is.null(multi_spec))) TRUE else if (is.null(multi_spec)) NULL else multi_spec
cox_result <- tryCatch(
.r4vn_surv_cox_models(
d_model, focal_map, map, adj_obj, mul_obj,
start = ".start" %in% names(d_model),
strata_name = if (".strata" %in% names(d_model)) ".strata" else NULL,
cluster_name = if (".cluster" %in% names(d_model)) ".cluster" else NULL,
frailty_name = if (".frailty" %in% names(d_model)) ".frailty" else NULL,
ci = ci, ties = ties, interaction_vars = interaction_vars
),
error = function(e) {
if (isTRUE(strict)) stop(e)
warning("Cox analysis was not available: ", conditionMessage(e), call. = FALSE)
list(crude = NULL, adjusted = NULL, multi = NULL,
interaction = NULL, diagnostics = NULL,
crude_fits = list(), adjusted_fits = list(), multi_fit = NULL)
}
)
if (!isTRUE(cox)) cox_result$crude <- NULL
}
ph_table <- NULL
if (isTRUE(ph)) {
ph_table <- tryCatch({
ph_fit <- cox_result$multi_fit
if (is.null(ph_fit)) {
focal_map <- map[match(focal_spec$variable, map$variable), , drop = FALSE]
ph_fit <- .r4vn_surv_fit_cox(
d_model, focal_map, ".start" %in% names(d_model),
if (".strata" %in% names(d_model)) ".strata" else NULL,
if (".cluster" %in% names(d_model)) ".cluster" else NULL,
if (".frailty" %in% names(d_model)) ".frailty" else NULL,
ties = ties
)
if (is.null(cox_result$multi_fit)) cox_result$multi_fit <- ph_fit
}
.r4vn_surv_ph(ph_fit)
}, error = function(e) {
if (isTRUE(strict)) stop(e)
warning("PH test unavailable: ", conditionMessage(e), call. = FALSE)
NULL
})
}
fg_result <- NULL
if (isTRUE(finegray)) {
if (!competing) stop("`finegray = TRUE` requires `compete`.", call. = FALSE)
fg_spec <- if (!is.null(multi_spec) && !isTRUE(multi_spec)) multi_spec else focal_spec
if (is.null(fg_spec)) stop("`finegray = TRUE` requires `vars`.", call. = FALSE)
fg_map <- map[match(fg_spec$variable, map$variable), , drop = FALSE]
fg_result <- tryCatch(
.r4vn_surv_finegray(d_model, fg_map, ".start" %in% names(d_model), ci, ties),
error = function(e) {
if (isTRUE(strict)) stop(e)
warning("Fine-Gray analysis was not available: ", conditionMessage(e), call. = FALSE)
NULL
}
)
}
# Subgroup final Cox models ---------------------------------------------
subgroup <- NULL
if (".superby" %in% names(d_model)) {
if (is.null(focal_spec) || !nrow(focal_spec)) stop("`superby` requires `vars` for subgroup Cox models.", call. = FALSE)
subgroup <- lapply(levels(d_model$.superby), function(g) {
z <- d_model[d_model$.superby == g, , drop = FALSE]
focal_map <- map[match(focal_spec$variable, map$variable), , drop = FALSE]
fit <- try(.r4vn_surv_fit_cox(z, focal_map, ".start" %in% names(z),
if (".strata" %in% names(z)) ".strata" else NULL,
if (".cluster" %in% names(z)) ".cluster" else NULL,
if (".frailty" %in% names(z)) ".frailty" else NULL,
ties = ties), silent = TRUE)
if (inherits(fit, "try-error")) return(list(group = g, fit = NULL, table = NULL))
list(group = g, fit = fit, table = .r4vn_surv_extract_cox(fit, focal_map, ci))
})
names(subgroup) <- levels(d_model$.superby)
}
# Compact analysis data kept for reusable plotting/risk tables.
analysis_data <- d_desc[, intersect(c(".time", ".event", ".start", ".by", ".id", ".status_ms"), names(d_desc)), drop = FALSE]
out <- list(
title = title,
overview = overview,
descriptive = descriptive,
followup = followup_table,
median = median_table,
lifetable = life_table,
at = at_table,
risk = if (isTRUE(risk) || isTRUE(rr) || isTRUE(rd) || !is.null(cuminc)) risk_table else NULL,
cuminc = if (isTRUE(risk) || isTRUE(rr) || isTRUE(rd) || !is.null(cuminc)) risk_table else NULL,
rate = rate_table,
risk_compare = risk_compare,
irr = irr_table,
logrank = logrank_table,
rmst = rmst_result,
cox = cox_result,
ph = ph_table,
finegray = fg_result,
subgroup = subgroup,
fit = fit,
curve = curve,
analysis_data = analysis_data,
metadata = list(time = time_name, event = event_name, by = by_name,
time_label = .r4vn_surv_label(data, time_name),
event_label = .r4vn_surv_label(data, event_name),
by_label = .r4vn_surv_label(data, by_name), id = id_name,
start = start_name, unit = unit, failure = status$failure, compete = compete,
vars = focal_spec, ci = ci, scale = scale, at = at,
cuminc = cuminc, competing = competing,
excluded = cc_desc$excluded, excluded_survival = cc_desc$excluded,
excluded_model_base = cc_model$excluded, missing = missing,
recurrent = recurrent, report = report,
lifetable = isTRUE(lifetable)),
call = call
)
class(out) <- "r4vn_surv"
if (isTRUE(ai)) out$ai_text <- .r4vn_surv_ai_text(out, effect_digit, p_digit)
out
}
.r4vn_tabsurv_by_info <- function(expr, data, env) {
if (.r4vn_expr_is_null(expr)) return(list(by = NULL, spec = NULL))
if (.r4vn_is_vars_spec_expr(expr, data, env)) {
spec <- .r4vn_by_spec(expr, data, env, allow_null = FALSE)
return(list(by = spec$by, spec = spec))
}
list(by = .r4vn_surv_name(expr, data, "by", TRUE), spec = NULL)
}
.r4vn_tabsurv_profile <- function(report, provided, values, data, time_name,
by_name = NULL, competing = FALSE,
recurrent = FALSE) {
n_vars <- if (is.null(values$vars)) 0L else if (inherits(values$vars, "r4vn_vars")) {
nrow(.r4vn_resolve_vars(values$vars, data = data, default_type = "auto", strict = TRUE))
} else if (is.character(values$vars)) length(unique(values$vars)) else 1L
has_vars <- n_vars > 0L
has_by <- !is.null(by_name)
n_groups <- if (has_by) {
length(unique(as.character(data[[by_name]][!is.na(data[[by_name]])])))
} else 0L
two_groups <- n_groups == 2L
legacy <- list(
followup = TRUE, km = TRUE, at = NULL, risk = FALSE, rate = FALSE,
logrank = TRUE, rr = FALSE, rd = FALSE, irr = FALSE,
cox = FALSE, adjusted = FALSE, multi = FALSE, finegray = FALSE,
rmst = FALSE, ph = FALSE, plot = FALSE
)
brief <- utils::modifyList(legacy, list(plot = TRUE))
auto <- list(
followup = TRUE, km = TRUE,
at = .r4vn_surv_auto_times(data[[time_name]]),
risk = TRUE, rate = "overall",
logrank = has_by && !competing,
rr = two_groups, rd = two_groups, irr = two_groups,
cox = has_vars, adjusted = FALSE, multi = n_vars > 1L,
finegray = competing && has_vars,
rmst = two_groups && !competing,
ph = has_vars, plot = TRUE
)
full <- utils::modifyList(auto, list(
adjusted = n_vars > 1L,
multi = has_vars,
rmst = has_by && !competing,
rate = "all"
))
defaults <- switch(report, custom = legacy, brief = brief, auto = auto, full = full)
out <- defaults
for (nm in names(out)) {
if (isTRUE(provided[[nm]]) && !is.null(values[[nm]])) out[[nm]] <- values[[nm]]
}
if (isTRUE(provided[["at"]])) out$at <- values$at
if (!isFALSE(recurrent)) {
out$cox <- TRUE
out$multi <- TRUE
for (nm in c("km", "risk", "logrank", "rr", "rd", "rmst", "plot", "finegray")) {
if (!isTRUE(provided[[nm]])) out[[nm]] <- FALSE
}
if (!isTRUE(provided[["at"]])) out$at <- NULL
}
if (competing) {
if (!isTRUE(provided[["logrank"]])) out$logrank <- FALSE
if (!isTRUE(provided[["rmst"]])) out$rmst <- FALSE
}
out
}
.r4vn_tabsurv_hierarchical_result <- function(results, labels, spec, outer,
data, call, title, missing,
ci, scale, report) {
n_total <- sum(vapply(results, function(z) z$overview$N[1L], numeric(1)), na.rm = TRUE)
ev_total <- sum(vapply(results, function(z) z$overview$Events[1L], numeric(1)), na.rm = TRUE)
cens_total <- sum(vapply(results, function(z) z$overview$Censored[1L], numeric(1)), na.rm = TRUE)
pt_total <- sum(vapply(results, function(z) z$overview$`Person-time`[1L], numeric(1)), na.rm = TRUE)
overview <- data.frame(
N = n_total, Events = ev_total, Censored = cens_total,
`Event percent` = if (n_total > 0L) 100 * ev_total / n_total else NA_real_,
`Person-time` = pt_total, check.names = FALSE
)
out <- list(
title = title, overview = overview,
hierarchical_results = results, strata_results = results,
strata_labels = labels,
hierarchical_by = list(
all = spec$all, strata = outer, by = spec$by,
all_labels = vapply(spec$all, function(nm) .r4vn_variable_label(data, nm), character(1)),
strata_labels = vapply(outer, function(nm) .r4vn_variable_label(data, nm), character(1)),
by_label = .r4vn_variable_label(data, spec$by)
),
metadata = list(
hierarchical = TRUE, by = spec$by, strata = outer, all_by = spec$all,
missing = missing, ci = ci, scale = scale, report = report
),
call = call
)
class(out) <- c("r4vn_surv_hierarchical", "r4vn_surv")
out
}
.r4vn_tabsurv_contract <- function(x, interpretation = FALSE) {
if (!inherits(x, "r4vn_surv")) return(x)
if (!is.null(x$hierarchical_results)) {
x$descriptive <- do.call(rbind, lapply(seq_along(x$hierarchical_results), function(i) {
z <- x$hierarchical_results[[i]]$descriptive
if (is.null(z) || !nrow(z)) return(NULL)
data.frame(Stratum = names(x$hierarchical_results)[i], z,
stringsAsFactors = FALSE, check.names = FALSE)
}))
x$estimates <- lapply(x$hierarchical_results, `[[`, "estimates")
x$tests <- lapply(x$hierarchical_results, `[[`, "tests")
x$diagnostics <- lapply(x$hierarchical_results, `[[`, "diagnostics")
x$models <- lapply(x$hierarchical_results, `[[`, "models")
if (isTRUE(interpretation)) {
x$interpretation <- do.call(rbind, lapply(seq_along(x$hierarchical_results), function(i) {
z <- x$hierarchical_results[[i]]$interpretation
if (is.null(z) || !nrow(z)) return(NULL)
data.frame(Stratum = names(x$hierarchical_results)[i], z,
stringsAsFactors = FALSE, check.names = FALSE)
}))
} else x$interpretation <- NULL
} else {
x$estimates <- list(
survival_at = x$at, life_table = x$lifetable,
cumulative_risk = x$risk, cumulative_incidence = x$cuminc,
incidence_rate = x$rate, rmst = x$rmst,
cox = x$cox[c("crude", "adjusted", "multi", "interaction")],
finegray = if (is.null(x$finegray)) NULL else x$finegray$table
)
x$tests <- list(
logrank = x$logrank, risk_comparison = x$risk_compare,
incidence_rate_ratio = x$irr
)
x$diagnostics <- list(
cox = if (is.null(x$cox)) NULL else x$cox$diagnostics,
proportional_hazards = x$ph
)
x$models <- list(
survival = x$fit,
crude_cox = if (is.null(x$cox)) list() else x$cox$crude_fits,
adjusted_cox = if (is.null(x$cox)) list() else x$cox$adjusted_fits,
multivariable_cox = if (is.null(x$cox)) NULL else x$cox$multi_fit,
finegray = if (is.null(x$finegray)) NULL else x$finegray$fit
)
if (isTRUE(interpretation)) {
x$interpretation <- .r4vn_surv_interpretation(x)
} else x$interpretation <- NULL
}
x$tables <- surv_tables(x)
x
}
.r4vn_tabsurv_plot <- function(x, plot = FALSE, plot_args = list(),
report = "custom", show = TRUE) {
if (!isTRUE(plot)) return(x)
gf <- get0("gsurv", mode = "function", inherits = TRUE)
if (is.null(gf)) {
warning("`plot = TRUE` requested but `gsurv()` is not available.", call. = FALSE)
return(x)
}
defaults <- if (identical(report, "custom")) list() else list(
ci = TRUE, pvalue = TRUE, risk_table = TRUE, median = FALSE
)
args <- utils::modifyList(defaults, plot_args)
if ("x" %in% names(args)) args$x <- NULL
if (is.null(args$show)) args$show <- isTRUE(show)
if (!isTRUE(show)) args$show <- FALSE
draw_one <- function(z, label = NULL) tryCatch({
args2 <- args
if (is.null(args2$legend) && !is.null(z$metadata$by_label)) {
args2$legend <- z$metadata$by_label
}
if (!is.null(label) && !is.null(args2$file)) {
ext <- tools::file_ext(args2$file)
stem <- if (nzchar(ext)) {
substr(args2$file, 1L, nchar(args2$file) - nchar(ext) - 1L)
} else args2$file
suffix <- gsub("[^A-Za-z0-9]+", "-", label)
suffix <- gsub("(^-+|-+$)", "", suffix)
args2$file <- paste0(stem, "-", suffix, if (nzchar(ext)) paste0(".", ext) else "")
}
do.call(gf, c(list(x = z), args2))
},
error = function(e) {
warning("Survival graph was not available: ", conditionMessage(e), call. = FALSE)
NULL
}
)
if (!is.null(x$hierarchical_results)) {
x$graphs <- Map(draw_one, x$hierarchical_results, names(x$hierarchical_results))
names(x$graphs) <- names(x$hierarchical_results)
x$plots <- x$graphs
} else {
x$graph <- draw_one(x)
x$plots <- list(survival = x$graph)
}
x
}
#' @export
print.r4vn_surv <- function(x, ...) {
if (!is.null(x$hierarchical_results)) {
if (!is.null(x$title) && nzchar(x$title)) cat(x$title, "\n", sep = "")
cat("R4VN hierarchical survival analysis\n")
cat(strrep("-", 58), "\n", sep = "")
cat("Hierarchy: ", paste(x$hierarchical_by$strata_labels %||%
x$hierarchical_by$strata, collapse = " > "),
" | Within-stratum group: ", x$hierarchical_by$by_label %||%
x$hierarchical_by$by, "\n", sep = "")
for (i in seq_along(x$hierarchical_results)) {
lab <- names(x$hierarchical_results)[i] %||% paste0("Stratum ", i)
cat("\n=== ", lab, " ===\n", sep = "")
print(x$hierarchical_results[[i]])
}
return(invisible(x))
}
if (!is.null(x$title) && nzchar(x$title)) cat(x$title, "\n", sep = "")
cat("R4VN survival analysis\n")
cat(strrep("-", 58), "\n", sep = "")
ov <- x$overview
cat(sprintf("N: %d | Events: %d | Censored: %d | Person-time: %.2f\n",
ov$N[1L], ov$Events[1L], ov$Censored[1L], ov$`Person-time`[1L]))
if (!is.null(x$descriptive) && nrow(x$descriptive)) {
cat("\nOutcome by ", x$metadata$by_label %||% "group", "\n", sep = "")
print(x$descriptive, row.names = FALSE)
}
if (isTRUE(x$metadata$missing) || x$metadata$excluded_survival > 0L || x$metadata$excluded_model_base > 0L) {
cat(sprintf("Excluded from survival summaries because time/event/group fields were incomplete: %d\n", x$metadata$excluded_survival))
if (x$metadata$excluded_model_base != x$metadata$excluded_survival) {
cat(sprintf("Excluded from the model base because structural model fields were incomplete: %d\n", x$metadata$excluded_model_base))
}
}
if (!is.null(x$followup) && nrow(x$followup)) {
cat("Median follow-up: ", .r4vn_surv_ci_text(x$followup$median[1], x$followup$lower[1], x$followup$upper[1], 2), "\n", sep = "")
}
if (!is.null(x$median) && nrow(x$median)) {
cat("\nMedian survival\n")
z <- x$median
z$estimate <- .r4vn_surv_ci_text(z$median, z$lower, z$upper, 2)
print(z[, c("group", "estimate"), drop = FALSE], row.names = FALSE)
}
if (!is.null(x$lifetable) && nrow(x$lifetable)) {
cat("\nLife table", if (isTRUE(x$metadata$competing))
" (Aalen-Johansen event history)" else " (Kaplan-Meier)", "\n", sep = "")
print(x$lifetable, row.names = FALSE)
}
if (!is.null(x$at) && nrow(x$at) && is.null(x$risk)) {
cat("\nSurvival probability at requested times\n")
z <- x$at
z$estimate <- .r4vn_surv_ci_text(z$surv, z$lower, z$upper, 2, percent = TRUE)
print(z[, intersect(c("group", "time", "n.risk", "estimate"), names(z)), drop = FALSE], row.names = FALSE)
}
if (!is.null(x$risk) && nrow(x$risk)) {
cat("\nCumulative risk", if (isTRUE(x$metadata$competing)) " (Aalen-Johansen CIF)" else " (1-KM)", "\n", sep = "")
z <- x$risk
z$estimate <- .r4vn_surv_ci_text(z$risk, z$risk_lower, z$risk_upper, 2, percent = TRUE)
print(z[, intersect(c("group", "time", "n.risk", "estimate"), names(z)), drop = FALSE], row.names = FALSE)
}
if (!is.null(x$rate) && nrow(x$rate)) {
cat("\nIncidence rate\n")
z <- x$rate
z$estimate <- .r4vn_surv_ci_text(z$rate, z$lower, z$upper, 2)
print(z[, c("group", "interval", "events", "person_time", "estimate"), drop = FALSE], row.names = FALSE)
}
if (!is.null(x$logrank)) cat("\nLog-rank: chi-square = ", .r4vn_surv_fmt(x$logrank$chisq[1], 2),
", df = ", x$logrank$df[1], ", p ", .r4vn_surv_fmt_p(x$logrank$p[1], 3), "\n", sep = "")
if (!is.null(x$risk_compare)) {
cat("\nRisk comparison\n"); print(x$risk_compare, row.names = FALSE)
}
if (!is.null(x$irr)) {
cat("\nIncidence rate ratio\n"); print(x$irr, row.names = FALSE)
}
print_model <- function(z, heading, effect = "HR") {
if (is.null(z) || !nrow(z)) return(invisible(NULL))
cat("\n", heading, "\n", sep = "")
zz <- z
zz[[paste0(effect, " (95% CI)")]] <- ifelse(zz$reference, "Ref.", .r4vn_surv_ci_text(zz$estimate, zz$lower, zz$upper, 2))
zz$p <- vapply(zz$p, .r4vn_surv_fmt_p, character(1), digits = 3)
print(zz[, intersect(c("variable", "variable_label", "level",
paste0(effect, " (95% CI)"), "p"), names(zz)), drop = FALSE],
row.names = FALSE)
}
print_model(x$cox$crude, "Crude Cox proportional hazards", "HR")
print_model(x$cox$adjusted, "Adjusted Cox proportional hazards", "HR")
print_model(x$cox$multi, "Multivariable Cox proportional hazards", "HR")
if (!is.null(x$cox$interaction) && nrow(x$cox$interaction)) {
cat("\nInteraction terms\n")
z <- x$cox$interaction
z$`HR (95% CI)` <- .r4vn_surv_ci_text(z$estimate, z$lower, z$upper, 2)
z$p <- vapply(z$p, .r4vn_surv_fmt_p, character(1), digits = 3)
print(z[, c("term", "HR (95% CI)", "p"), drop = FALSE], row.names = FALSE)
}
if (!is.null(x$cox$diagnostics) && nrow(x$cox$diagnostics)) {
z <- x$cox$diagnostics[1L, ]
cat(sprintf("\nModel N = %d; events = %d; concordance = %.3f; LR p %s\n",
as.integer(z$n), as.integer(z$events), z$concordance, .r4vn_surv_fmt_p(z$LR_p, 3)))
}
if (!is.null(x$finegray)) print_model(x$finegray$table, "Fine-Gray subdistribution hazards", "SHR")
if (!is.null(x$ph)) {
cat("\nProportional-hazards test\n")
z <- x$ph
if ("p" %in% names(z)) z$p <- vapply(z$p, .r4vn_surv_fmt_p, character(1), digits = 3)
print(z, row.names = FALSE)
}
if (!is.null(x$rmst)) {
cat("\nRestricted mean survival time\n")
print(x$rmst$table, row.names = FALSE)
if (!is.null(x$rmst$difference)) print(x$rmst$difference, row.names = FALSE)
}
if (!is.null(x$interpretation) && nrow(x$interpretation)) {
cat("\nInterpretation\n")
for (i in seq_len(nrow(x$interpretation))) {
cat("- ", x$interpretation$section[i], ": ",
x$interpretation$interpretation[i], "\n", sep = "")
}
}
invisible(x)
}
#' Direct Cox Proportional Hazards Model
#'
#' A compact R4VN wrapper around `tabsurv()` for a final multivariable Cox model.
#'
#' @inheritParams tabsurv
#' @param ph Logical; test the proportional-hazards assumption. Default `FALSE` in `cox()`. Setting `diagnosis = TRUE` also requests this diagnostic.
#' @param diagnosis Logical; if `TRUE`, append Cox-model diagnostics including concordance, proportional-hazards testing, and residual summaries. Default `FALSE`.
#' @return An object of class `r4vn_surv`.
#' @examples
#' if (requireNamespace("survival", quietly = TRUE)) {
#' d <- data.frame(
#' time = c(5, 8, 10, 12, 15, 18, 20, 22, 25, 30),
#' event = c(1, 0, 1, 1, 0, 1, 0, 1, 1, 0),
#' age = c(40, 45, 50, 55, 60, 48, 52, 63, 58, 67),
#' sex = factor(rep(c("Female", "Male"), 5))
#' )
#' cox(time, event, vars = vars(c.age, i.sex), data = d, show = FALSE)
#' cox(time, event, vars = vars(c.age, i.sex), data = d,
#' diagnosis = TRUE, show = FALSE)
#' cox(time, event, vars = vars(c.age), data = d, ph = TRUE, show = FALSE)
#' }
#' @export
cox <- function(time, event, vars, data = NULL, failure = NULL, id = NULL,
start = NULL, strata = NULL, cluster = NULL, frailty = NULL,
ph = FALSE, ci = .95, ties = c("efron", "breslow", "exact"),
diagnosis = FALSE, show = TRUE, console = FALSE) {
call <- match.call()
# Resolve the default choice before forwarding to tabsurv().
ties <- match.arg(ties)
d <- .r4vn_surv_data(data)
tn <- .r4vn_surv_name(substitute(time), d, "time")
en <- .r4vn_surv_name(substitute(event), d, "event")
idn <- .r4vn_surv_name(substitute(id), d, "id", TRUE)
stn <- .r4vn_surv_name(substitute(start), d, "start", TRUE)
srn <- .r4vn_surv_name(substitute(strata), d, "strata", TRUE)
cln <- .r4vn_surv_name(substitute(cluster), d, "cluster", TRUE)
frn <- .r4vn_surv_name(substitute(frailty), d, "frailty", TRUE)
# tabsurv() uses NSE for time/event/id/start/strata/cluster/frailty. A
# direct call such as `id = idn` would therefore expose the symbol `idn` to
# tabsurv(), even when its value is NULL, and tabsurv() would look for a
# column literally named "idn". Build the call from evaluated values so
# character variable names and literal NULLs arrive intact.
args <- list(
time = tn, event = en, vars = vars, data = d, failure = failure,
id = idn, start = stn, strata = srn, cluster = cln, frailty = frn,
km = FALSE, cox = TRUE, multi = TRUE, ph = isTRUE(ph) || isTRUE(diagnosis), ci = ci,
ties = ties, report = "custom", plot = FALSE,
interpretation = FALSE, show = FALSE, console = FALSE
)
out <- do.call(tabsurv, args)
if (!isTRUE(diagnosis) && is.list(out$cox)) out$cox$diagnostics <- NULL
if (isTRUE(diagnosis) && !is.null(out$cox$multi_fit)) {
out$model_diagnostics <- .r4vn_model_diagnosis(out$cox$multi_fit, kind = "cox")
} else if (isTRUE(diagnosis) && !is.null(out$raw$model)) {
out$model_diagnostics <- .r4vn_model_diagnosis(out$raw$model, kind = "cox")
}
out$call <- call
.r4vn_show(out, show = show, console = console)
}
# ============================================================================
# Viewer renderers for survival commands
# ============================================================================
.r4vn_surv_model_display <- function(z, effect = "HR", digit = 2, p_digit = 3) {
if (is.null(z) || !nrow(z)) return(NULL)
zz <- z
zz[[paste0(effect, " (95% CI)")]] <- ifelse(zz$reference, "Ref.",
.r4vn_surv_ci_text(zz$estimate, zz$lower, zz$upper, digit))
zz$p <- vapply(zz$p, .r4vn_surv_fmt_p, character(1), digits = p_digit)
zz[, intersect(c("variable", "variable_label", "level", "term",
paste0(effect, " (95% CI)"), "p"), names(zz)), drop = FALSE]
}
.r4vn_surv_viewer <- function(x, subtitle = "Survival analysis") {
blocks <- character()
ov <- x$overview
if (!is.null(ov) && nrow(ov)) {
overview <- data.frame(Statistic = names(ov), Value = as.character(ov[1, ]), stringsAsFactors = FALSE)
blocks <- c(blocks, paste0('<section class="r4vn-section"><h2>Overview</h2>',
.r4vn_view_key_values(overview), '</section>'))
}
if (!is.null(x$descriptive) && nrow(x$descriptive)) {
blocks <- c(blocks, .r4vn_view_section(
paste0("Outcome by ", x$metadata$by_label %||% "group"), x$descriptive
))
}
if (!is.null(x$followup) && nrow(x$followup)) {
z <- x$followup; z$`Median follow-up (95% CI)` <- .r4vn_surv_ci_text(z$median, z$lower, z$upper, 2)
blocks <- c(blocks, .r4vn_view_section("Follow-up", z[, intersect(c("group", "Median follow-up (95% CI)"), names(z)), drop = FALSE]))
}
if (!is.null(x$median) && nrow(x$median)) {
z <- x$median; z$`Median survival (95% CI)` <- .r4vn_surv_ci_text(z$median, z$lower, z$upper, 2)
blocks <- c(blocks, .r4vn_view_section("Median survival", z[, intersect(c("group", "Median survival (95% CI)"), names(z)), drop = FALSE]))
}
if (!is.null(x$lifetable) && nrow(x$lifetable)) {
ttl <- if (isTRUE(x$metadata$competing)) {
"Life table (Aalen-Johansen event history)"
} else "Life table (Kaplan-Meier)"
blocks <- c(blocks, .r4vn_view_section(ttl, x$lifetable))
}
if (!is.null(x$at) && nrow(x$at) && is.null(x$risk)) {
z <- x$at; z$`Survival (95% CI)` <- .r4vn_surv_ci_text(z$surv, z$lower, z$upper, 2, percent = TRUE)
blocks <- c(blocks, .r4vn_view_section("Survival probability at requested times", z[, intersect(c("group", "time", "n.risk", "Survival (95% CI)"), names(z)), drop = FALSE]))
}
if (!is.null(x$risk) && nrow(x$risk)) {
z <- x$risk; z$`Risk (95% CI)` <- .r4vn_surv_ci_text(z$risk, z$risk_lower, z$risk_upper, 2, percent = TRUE)
ttl <- if (isTRUE(x$metadata$competing)) "Cumulative incidence (Aalen-Johansen CIF)" else "Cumulative risk (1 - KM)"
blocks <- c(blocks, .r4vn_view_section(ttl, z[, intersect(c("group", "time", "n.risk", "Risk (95% CI)"), names(z)), drop = FALSE]))
}
if (!is.null(x$rate) && nrow(x$rate)) {
z <- x$rate; z$`Rate (95% CI)` <- .r4vn_surv_ci_text(z$rate, z$lower, z$upper, 2)
blocks <- c(blocks, .r4vn_view_section("Incidence rate", z[, intersect(c("group", "interval", "events", "person_time", "Rate (95% CI)"), names(z)), drop = FALSE]))
}
if (!is.null(x$logrank) && nrow(x$logrank)) {
z <- data.frame(
`Chi-square` = .r4vn_surv_fmt(x$logrank$chisq, 2),
df = .r4vn_surv_fmt(x$logrank$df, 0),
p = vapply(x$logrank$p, .r4vn_surv_fmt_p, character(1), digits = 3),
stringsAsFactors = FALSE, check.names = FALSE
)
blocks <- c(blocks, .r4vn_view_section("Log-rank test", z))
}
if (!is.null(x$risk_compare) && nrow(x$risk_compare)) {
z <- x$risk_compare
out <- z[, intersect(c("time", "reference", "comparison"), names(z)), drop = FALSE]
if ("rr" %in% names(z)) {
out$`RR (95% CI)` <- .r4vn_surv_ci_text(z$rr, z$rr_lower, z$rr_upper, 2)
out$`RR p` <- vapply(z$rr_p, .r4vn_surv_fmt_p, character(1), digits = 3)
}
if ("rd" %in% names(z)) {
out$`RD (95% CI)` <- .r4vn_surv_ci_text(z$rd, z$rd_lower, z$rd_upper, 3)
out$`RD p` <- vapply(z$rd_p, .r4vn_surv_fmt_p, character(1), digits = 3)
}
blocks <- c(blocks, .r4vn_view_section("Risk comparison", out))
}
if (!is.null(x$irr) && nrow(x$irr)) {
z <- x$irr
out <- z[, intersect(c("interval", "reference", "comparison"), names(z)), drop = FALSE]
out$`IRR (95% CI)` <- .r4vn_surv_ci_text(z$irr, z$lower, z$upper, 2)
out$p <- vapply(z$p, .r4vn_surv_fmt_p, character(1), digits = 3)
blocks <- c(blocks, .r4vn_view_section("Incidence rate ratio", out))
}
if (!is.null(x$cox$crude)) blocks <- c(blocks, .r4vn_view_section("Crude Cox proportional hazards", .r4vn_surv_model_display(x$cox$crude, "HR")))
if (!is.null(x$cox$adjusted)) blocks <- c(blocks, .r4vn_view_section("Adjusted Cox proportional hazards", .r4vn_surv_model_display(x$cox$adjusted, "HR")))
if (!is.null(x$cox$multi)) blocks <- c(blocks, .r4vn_view_section("Multivariable Cox proportional hazards", .r4vn_surv_model_display(x$cox$multi, "HR")))
if (!is.null(x$cox$interaction) && nrow(x$cox$interaction)) {
z <- x$cox$interaction; z$`HR (95% CI)` <- .r4vn_surv_ci_text(z$estimate, z$lower, z$upper, 2); z$p <- vapply(z$p, .r4vn_surv_fmt_p, character(1), digits = 3)
blocks <- c(blocks, .r4vn_view_section("Interaction terms", z[, intersect(c("term", "HR (95% CI)", "p"), names(z)), drop = FALSE]))
}
if (!is.null(x$cox$diagnostics) && nrow(x$cox$diagnostics)) {
z <- x$cox$diagnostics
pcols <- grep("_p$|^p$", names(z), ignore.case = TRUE, value = TRUE)
for (nm in pcols) z[[nm]] <- vapply(z[[nm]], .r4vn_surv_fmt_p, character(1), digits = 3)
for (nm in setdiff(names(z), pcols)) {
if (is.numeric(z[[nm]])) z[[nm]] <- if (nm %in% c("n", "events", "LR_df", "Wald_df", "Score_df")) .r4vn_surv_fmt(z[[nm]], 0) else .r4vn_surv_fmt(z[[nm]], 3)
}
blocks <- c(blocks, .r4vn_view_section("Cox model diagnostics", z))
}
if (!is.null(x$finegray$table)) blocks <- c(blocks, .r4vn_view_section("Fine-Gray subdistribution hazards", .r4vn_surv_model_display(x$finegray$table, "SHR")))
if (!is.null(x$model_diagnostics) && length(x$model_diagnostics)) {
for (nm in names(x$model_diagnostics)) {
z <- x$model_diagnostics[[nm]]
if (is.data.frame(z) && nrow(z)) blocks <- c(blocks, .r4vn_view_section(nm, z))
}
}
if (!is.null(x$ph) && nrow(x$ph)) {
z <- x$ph
for (nm in names(z)) {
if (tolower(nm) == "p") z[[nm]] <- vapply(z[[nm]], .r4vn_surv_fmt_p, character(1), digits = 3)
else if (is.numeric(z[[nm]])) z[[nm]] <- .r4vn_surv_fmt(z[[nm]], 3)
}
blocks <- c(blocks, .r4vn_view_section("Proportional-hazards test", z))
}
if (!is.null(x$rmst$table) && nrow(x$rmst$table)) {
z <- x$rmst$table
out <- z[, intersect(c("group", "tau"), names(z)), drop = FALSE]
if ("tau" %in% names(out)) out$tau <- .r4vn_surv_fmt(out$tau, 2)
out$`RMST (95% CI)` <- .r4vn_surv_ci_text(z$rmst, z$lower, z$upper, 2)
out$SE <- .r4vn_surv_fmt(z$se, 2)
blocks <- c(blocks, .r4vn_view_section("Restricted mean survival time", out))
}
if (!is.null(x$rmst$difference) && nrow(x$rmst$difference)) {
z <- x$rmst$difference
out <- z[, intersect(c("reference", "comparison"), names(z)), drop = FALSE]
out$`Difference (95% CI)` <- .r4vn_surv_ci_text(z$difference, z$lower, z$upper, 2)
out$p <- vapply(z$p, .r4vn_surv_fmt_p, character(1), digits = 3)
blocks <- c(blocks, .r4vn_view_section("RMST difference", out))
}
if (!is.null(x$interpretation) && nrow(x$interpretation)) {
blocks <- c(blocks, .r4vn_view_section("Interpretation", x$interpretation))
}
notes <- character()
if (isTRUE(x$metadata$missing) || x$metadata$excluded_survival > 0L)
notes <- c(notes, paste0("Excluded from survival summaries because required fields were incomplete: ", x$metadata$excluded_survival, "."))
if (!is.null(x$metadata$excluded_model_base) && x$metadata$excluded_model_base != x$metadata$excluded_survival)
notes <- c(notes, paste0("Excluded from the structural model base: ", x$metadata$excluded_model_base, "."))
.r4vn_view_document(x$title %||% "R4VN survival analysis", paste0(blocks, collapse = ""), notes = notes,
subtitle = subtitle, prefix = "r4vn-tabsurv-")
}
.r4vn_surv_hierarchical_viewer <- function(x) {
results <- x$hierarchical_results
labels <- names(results) %||% x$strata_labels %||% paste0("Stratum ", seq_along(results))
summary_rows <- lapply(seq_along(results), function(i) {
ov <- results[[i]]$overview
data.frame(
Stratum = labels[i], N = ov$N[1L], Events = ov$Events[1L],
Censored = ov$Censored[1L], `Event %` = .r4vn_surv_fmt(ov$`Event percent`[1L], 1),
`Person-time` = .r4vn_surv_fmt(ov$`Person-time`[1L], 2),
stringsAsFactors = FALSE, check.names = FALSE
)
})
blocks <- .r4vn_view_section("Hierarchy overview", do.call(rbind, summary_rows))
for (i in seq_along(results)) {
z <- results[[i]]
blocks <- paste0(
blocks, '<section class="r4vn-section"><h2>',
.r4vn_view_escape(labels[i]), '</h2></section>'
)
if (!is.null(z$descriptive) && nrow(z$descriptive)) {
blocks <- paste0(blocks, .r4vn_view_section("Outcome by group", z$descriptive))
}
if (!is.null(z$median) && nrow(z$median)) {
q <- z$median
q$`Median survival (95% CI)` <- .r4vn_surv_ci_text(q$median, q$lower, q$upper, 2)
blocks <- paste0(blocks, .r4vn_view_section(
"Median survival", q[, intersect(c("group", "Median survival (95% CI)"), names(q)), drop = FALSE]
))
}
if (!is.null(z$lifetable) && nrow(z$lifetable)) {
blocks <- paste0(blocks, .r4vn_view_section(
if (isTRUE(z$metadata$competing)) {
"Life table (Aalen-Johansen event history)"
} else "Life table (Kaplan-Meier)",
z$lifetable
))
}
if (!is.null(z$risk) && nrow(z$risk)) {
q <- z$risk
q$`Risk (95% CI)` <- .r4vn_surv_ci_text(q$risk, q$risk_lower, q$risk_upper, 2, percent = TRUE)
blocks <- paste0(blocks, .r4vn_view_section(
if (isTRUE(z$metadata$competing)) "Cumulative incidence" else "Cumulative risk",
q[, intersect(c("group", "time", "n.risk", "Risk (95% CI)"), names(q)), drop = FALSE]
))
}
if (!is.null(z$rate) && nrow(z$rate)) {
q <- z$rate
q$`Rate (95% CI)` <- .r4vn_surv_ci_text(q$rate, q$lower, q$upper, 2)
blocks <- paste0(blocks, .r4vn_view_section(
"Incidence rate", q[, intersect(c("group", "interval", "events", "person_time", "Rate (95% CI)"), names(q)), drop = FALSE]
))
}
if (!is.null(z$logrank) && nrow(z$logrank)) {
q <- data.frame(
`Chi-square` = .r4vn_surv_fmt(z$logrank$chisq, 2),
df = .r4vn_surv_fmt(z$logrank$df, 0),
p = vapply(z$logrank$p, .r4vn_surv_fmt_p, character(1), digits = 3),
stringsAsFactors = FALSE, check.names = FALSE
)
blocks <- paste0(blocks, .r4vn_view_section("Log-rank test", q))
}
model <- if (!is.null(z$cox$multi)) z$cox$multi else z$cox$crude
if (!is.null(model) && nrow(model)) {
blocks <- paste0(blocks, .r4vn_view_section(
if (!is.null(z$cox$multi)) "Multivariable Cox proportional hazards" else "Crude Cox proportional hazards",
.r4vn_surv_model_display(model, "HR")
))
}
if (!is.null(z$interpretation) && nrow(z$interpretation)) {
blocks <- paste0(blocks, .r4vn_view_section("Interpretation", z$interpretation))
}
}
notes <- paste0(
"Hierarchical by: ", paste(x$hierarchical_by$strata_labels %||%
x$hierarchical_by$strata, collapse = " > "),
"; innermost survival group: ", x$hierarchical_by$by_label %||%
x$hierarchical_by$by, "."
)
.r4vn_view_document(
x$title %||% "R4VN hierarchical survival analysis", blocks,
notes = notes, subtitle = "Hierarchical survival analysis",
prefix = "r4vn-tabsurv-hierarchical-"
)
}
.r4vn_viewer_tabsurv <- function(x) {
if (!is.null(x$hierarchical_results)) return(.r4vn_surv_hierarchical_viewer(x))
.r4vn_surv_viewer(x, "Survival, incidence and time-to-event analysis")
}
.r4vn_viewer_cox <- function(x) .r4vn_surv_viewer(x, "Cox proportional hazards model")
Try the R4VN package in your browser
Any scripts or data that you put into this service are public.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.