R/surv-utils.R

In R4VN: Health Data Analysis and Publication-Ready Reporting

# R4VN survival utilities -------------------------------------------------
# Internal helpers for tabsurv(), cox(), and gsurv().

.r4vn_surv_require <- function() {
  if (!requireNamespace("survival", quietly = TRUE)) {
    stop("Package `survival` is required. Install it with install.packages('survival').", call. = FALSE)
  }
  invisible(TRUE)
}

.r4vn_surv_data <- function(data = NULL) {
  if (is.null(data)) {
    f <- get0(".r4vn_get_active", mode = "function", inherits = TRUE)
    if (!is.null(f)) return(f())
    f <- get0(".r4vn_resolve_analysis_data", mode = "function", inherits = TRUE)
    if (!is.null(f)) return(f(NULL))
    stop("No active data frame. Use `usedf(data)` or supply `data =`.", call. = FALSE)
  }
  if (!is.data.frame(data)) stop("`data` must be a data frame.", call. = FALSE)
  data
}

.r4vn_surv_deparse <- function(x) paste(deparse(x, width.cutoff = 500L), collapse = "")

.r4vn_surv_name <- function(expr, data, arg, allow_null = FALSE) {
  if (identical(expr, quote(NULL)) || identical(expr, quote(expr = ))) {
    if (allow_null) return(NULL)
    stop(sprintf("`%s` is required.", arg), call. = FALSE)
  }
  if (is.symbol(expr)) {
    nm <- as.character(expr)
    if (nm %in% names(data)) return(nm)
  }
  val <- try(eval(expr, envir = parent.frame(2L)), silent = TRUE)
  if (!inherits(val, "try-error") && is.character(val) && length(val) == 1L && val %in% names(data)) return(val)
  txt <- .r4vn_surv_deparse(expr)
  stop(sprintf("`%s` must identify one variable in `data`; `%s` was not found.", arg, txt), call. = FALSE)
}

.r4vn_surv_label <- function(data, name) {
  if (is.null(name)) return(NULL)
  lab <- attr(data[[name]], "label", exact = TRUE)
  if (!is.null(lab) && length(lab) && nzchar(as.character(lab)[1L])) as.character(lab)[1L] else name
}

.r4vn_surv_unit_label <- function(unit, plural = TRUE) {
  if (is.null(unit) || !length(unit) || is.na(unit[1L]) || !nzchar(as.character(unit)[1L])) return("time units")
  z <- tolower(as.character(unit)[1L])
  dict <- c(day = "day", days = "day", month = "month", months = "month", year = "year", years = "year",
            week = "week", weeks = "week", hour = "hour", hours = "hour")
  ans <- if (z %in% names(dict)) unname(dict[z]) else as.character(unit)[1L]
  if (plural) paste0(ans, "s") else ans
}

.r4vn_surv_fmt_p <- function(p, digits = 3L) {
  if (length(p) != 1L || is.na(p) || !is.finite(p)) return(NA_character_)
  cut <- 10^(-digits)
  if (p < cut) return(paste0("<", formatC(cut, format = "f", digits = digits)))
  formatC(p, format = "f", digits = digits)
}

.r4vn_surv_fmt <- function(x, digits = 2L) {
  ifelse(is.na(x) | !is.finite(x), NA_character_, formatC(x, format = "f", digits = digits))
}

.r4vn_surv_ci_text <- function(est, lo, hi, digits = 2L, percent = FALSE) {
  if (percent) {
    est <- est * 100; lo <- lo * 100; hi <- hi * 100
  }
  paste0(.r4vn_surv_fmt(est, digits), " (", .r4vn_surv_fmt(lo, digits), "-", .r4vn_surv_fmt(hi, digits), ")")
}

.r4vn_surv_as_spec <- function(x, data, arg = "vars", allow_null = TRUE) {
  if (is.null(x)) {
    if (allow_null) return(NULL)
    stop(sprintf("`%s` is required.", arg), call. = FALSE)
  }
  if (inherits(x, "r4vn_vars")) {
    resolver <- get0(".r4vn_resolve_vars", mode = "function", inherits = TRUE)
    out <- if (is.null(resolver)) {
      as.data.frame(x, stringsAsFactors = FALSE)
    } else {
      as.data.frame(
        resolver(x, data = data, default_type = "auto", strict = TRUE),
        stringsAsFactors = FALSE
      )
    }
    need <- c("variable", "type", "reference_index")
    if (!all(need %in% names(out))) stop(sprintf("`%s` is not a compatible `vars()` object.", arg), call. = FALSE)
    if (!"specification" %in% names(out)) out$specification <- out$variable
  } else if (is.character(x)) {
    if (!length(x)) return(NULL)
    typ <- vapply(x, function(nm) if (nm %in% names(data) && is.numeric(data[[nm]])) "mean" else "categorical", character(1))
    out <- data.frame(variable = x, type = typ, specification = x,
                      reference_index = ifelse(typ == "categorical", 1L, NA_integer_), stringsAsFactors = FALSE)
  } else {
    stop(sprintf("`%s` must be created using `vars()` or be a character vector of variable names.", arg), call. = FALSE)
  }
  miss <- setdiff(out$variable, names(data))
  if (length(miss)) stop(sprintf("Variables not found in `data`: %s.", paste(miss, collapse = ", ")), call. = FALSE)
  out <- out[!duplicated(out$variable), , drop = FALSE]
  rownames(out) <- NULL
  out
}

.r4vn_surv_merge_specs <- function(...) {
  xs <- Filter(Negate(is.null), list(...))
  if (!length(xs)) return(NULL)
  z <- do.call(rbind, xs)
  z <- z[!duplicated(z$variable), , drop = FALSE]
  rownames(z) <- NULL
  z
}

.r4vn_surv_apply_reference <- function(x, spec_row) {
  if (!identical(spec_row$type, "categorical")) {
    if (!is.numeric(x)) stop(sprintf("`%s` must be numeric because it was specified with c./q./f.", spec_row$variable), call. = FALSE)
    return(as.numeric(x))
  }
  if (is.factor(x)) {
    f <- droplevels(x)
  } else {
    lev <- unique(as.character(x[!is.na(x)]))
    f <- factor(as.character(x), levels = lev)
  }
  idx <- spec_row$reference_index
  if (is.na(idx)) idx <- 1L
  if (idx < 1L || idx > nlevels(f)) {
    stop(sprintf("Reference b%d.%s is invalid because `%s` has %d observed levels.", idx, spec_row$variable, spec_row$variable, nlevels(f)), call. = FALSE)
  }
  stats::relevel(f, ref = levels(f)[idx])
}

.r4vn_surv_internal_data <- function(data, spec = NULL) {
  d <- data.frame(row.names = seq_len(nrow(data)))
  map <- NULL
  if (!is.null(spec) && nrow(spec)) {
    map <- spec
    map$internal <- paste0(".x", seq_len(nrow(map)))
    map$variable_label <- vapply(
      map$variable, function(nm) .r4vn_surv_label(data, nm), character(1)
    )
    for (i in seq_len(nrow(map))) d[[map$internal[i]]] <- .r4vn_surv_apply_reference(data[[map$variable[i]]], map[i, , drop = FALSE])
  }
  list(data = d, map = map)
}

.r4vn_surv_default_failure <- function(x) {
  z <- x[!is.na(x)]
  if (!length(z)) stop("The event variable contains no non-missing values.", call. = FALSE)
  if (is.logical(z)) return(TRUE)
  if (is.factor(x)) {
    lev <- levels(droplevels(x))
    if (length(lev) != 2L) stop("`failure` must be supplied when the event variable has more than two levels.", call. = FALSE)
    return(lev[2L])
  }
  u <- unique(z)
  if (length(u) != 2L) stop("`failure` must be supplied when the event variable does not have exactly two observed values.", call. = FALSE)
  if (is.numeric(u)) {
    if (all(sort(u) == c(0, 1))) return(1)
    return(max(u))
  }
  as.character(u[2L])
}

.r4vn_surv_equal <- function(x, value) {
  if (is.factor(x)) x <- as.character(x)
  if (is.factor(value)) value <- as.character(value)
  if (is.character(x) || is.character(value)) return(as.character(x) %in% as.character(value))
  x %in% value
}

.r4vn_surv_event_status <- function(x, failure = NULL, compete = NULL) {
  if (is.null(failure)) failure <- .r4vn_surv_default_failure(x)
  fail <- .r4vn_surv_equal(x, failure)
  comp <- if (is.null(compete)) rep(FALSE, length(x)) else .r4vn_surv_equal(x, compete)
  fail[is.na(x)] <- NA
  comp[is.na(x)] <- NA
  if (any(fail & comp, na.rm = TRUE)) stop("`failure` and `compete` overlap.", call. = FALSE)

  out01 <- ifelse(is.na(x), NA_integer_, as.integer(fail))
  ms <- rep("censor", length(x))
  ms[fail %in% TRUE] <- "failure"
  if (!is.null(compete)) {
    raw <- if (is.factor(x)) as.character(x) else as.character(x)
    for (v in as.character(compete)) ms[!is.na(raw) & raw == v] <- paste0("compete:", v)
  }
  ms[is.na(x)] <- NA_character_
  lev <- c("censor", "failure", sort(unique(ms[grepl("^compete:", ms)])))
  status_ms <- factor(ms, levels = lev)
  list(event = out01, failure = failure, compete = compete, status_ms = status_ms)
}

.r4vn_surv_validate_time <- function(time, event, start = NULL) {
  if (!is.numeric(time)) stop("`time` must be numeric.", call. = FALSE)
  if (!is.null(start) && !is.numeric(start)) stop("`start` must be numeric.", call. = FALSE)
  if (any(time < 0, na.rm = TRUE)) stop("`time` cannot contain negative values.", call. = FALSE)
  if (!is.null(start)) {
    if (any(start < 0, na.rm = TRUE)) stop("`start` cannot contain negative values.", call. = FALSE)
    bad <- !is.na(start) & !is.na(time) & start >= time
    if (any(bad)) stop("Every non-missing `start` value must be smaller than `time` (stop time).", call. = FALSE)
  }
  if (!any(event == 1L, na.rm = TRUE)) warning("No events of interest were observed.", call. = FALSE)
  invisible(TRUE)
}

.r4vn_surv_complete <- function(d, columns) {
  columns <- unique(columns[!is.na(columns) & nzchar(columns)])
  keep <- stats::complete.cases(d[, columns, drop = FALSE])
  list(data = d[keep, , drop = FALSE], keep = keep, excluded = sum(!keep))
}

.r4vn_surv_formula <- function(start = FALSE, rhs = "1") {
  lhs <- if (isTRUE(start)) "survival::Surv(.start, .time, .event)" else "survival::Surv(.time, .event)"
  stats::as.formula(paste(lhs, "~", rhs), env = parent.frame())
}

.r4vn_surv_ms_formula <- function(start = FALSE, rhs = "1") {
  lhs <- if (isTRUE(start)) "survival::Surv(.start, .time, .status_ms)" else "survival::Surv(.time, .status_ms)"
  stats::as.formula(paste(lhs, "~", rhs), env = parent.frame())
}

.r4vn_surv_clean_strata <- function(x) {
  if (is.null(x)) return(rep("All", 0L))
  z <- as.character(x)
  z <- sub("^\\.by=", "", z)
  z <- sub("^by=", "", z)
  z
}

.r4vn_surv_curve_km <- function(fit) {
  if (is.null(fit$strata)) {
    idxs <- list(All = seq_along(fit$time))
  } else {
    ends <- cumsum(as.integer(fit$strata)); starts <- c(1L, head(ends, -1L) + 1L)
    nms <- .r4vn_surv_clean_strata(names(fit$strata))
    idxs <- Map(seq.int, starts, ends); names(idxs) <- nms
  }
  out <- lapply(names(idxs), function(g) {
    ii <- idxs[[g]]
    se <- if (is.null(fit$std.err)) rep(NA_real_, length(ii)) else fit$std.err[ii]
    data.frame(group = g, time = c(0, fit$time[ii]), estimate = c(1, fit$surv[ii]),
               std.err = c(0, se),
               lower = c(1, fit$lower[ii]), upper = c(1, fit$upper[ii]),
               n.risk = c(if (length(ii)) fit$n.risk[ii][1L] else NA_real_, fit$n.risk[ii]),
               n.event = c(0, fit$n.event[ii]), n.censor = c(0, fit$n.censor[ii]),
               kind = "survival", stringsAsFactors = FALSE)
  })
  do.call(rbind, out)
}

.r4vn_surv_lifetable <- function(curve, competing = FALSE) {
  groups <- unique(curve$group)
  out <- lapply(groups, function(g) {
    z <- curve[curve$group == g, , drop = FALSE]
    # The first row of every stored curve is the artificial time-zero row.
    if (nrow(z) <= 1L) return(NULL)
    z <- z[-1L, , drop = FALSE]
    base <- z[, intersect(c("group", "time", "n.risk", "n.event", "n.censor"),
                          names(z)), drop = FALSE]
    base$conditional_survival <- ifelse(
      is.finite(z$n.risk) & z$n.risk > 0,
      1 - z$n.event / z$n.risk,
      NA_real_
    )
    if (isTRUE(competing)) {
      base$cumulative_incidence <- z$estimate
    } else {
      base$survival <- z$estimate
      base$cumulative_risk <- 1 - z$estimate
    }
    base$std.error <- if ("std.err" %in% names(z)) z$std.err else NA_real_
    base$lower <- z$lower
    base$upper <- z$upper
    base
  })
  out <- Filter(Negate(is.null), out)
  if (!length(out)) return(data.frame())
  ans <- do.call(rbind, out)
  rownames(ans) <- NULL
  ans
}

.r4vn_surv_state_col <- function(fit, state = "failure") {
  cn <- colnames(fit$pstate)
  if (is.null(cn)) stop("Could not identify states in the Aalen-Johansen fit.", call. = FALSE)
  hit <- which(cn == state)
  if (!length(hit)) hit <- grep(paste0("(^|=)", state, "$"), cn)
  if (!length(hit)) hit <- grep(state, cn, fixed = TRUE)
  if (!length(hit)) stop(sprintf("State `%s` was not found in the Aalen-Johansen fit.", state), call. = FALSE)
  hit[1L]
}

.r4vn_surv_curve_cif <- function(fit, state = "failure") {
  j <- .r4vn_surv_state_col(fit, state)
  if (is.null(fit$strata)) {
    idxs <- list(All = seq_along(fit$time))
  } else {
    ends <- cumsum(as.integer(fit$strata)); starts <- c(1L, head(ends, -1L) + 1L)
    nms <- .r4vn_surv_clean_strata(names(fit$strata))
    idxs <- Map(seq.int, starts, ends); names(idxs) <- nms
  }
  out <- lapply(names(idxs), function(g) {
    ii <- idxs[[g]]
    lo <- if (!is.null(fit$lower)) fit$lower[ii, j] else NA_real_
    hi <- if (!is.null(fit$upper)) fit$upper[ii, j] else NA_real_
    se <- if (!is.null(fit$std.err)) {
      if (is.matrix(fit$std.err)) fit$std.err[ii, j] else fit$std.err[ii]
    } else rep(NA_real_, length(ii))
    data.frame(group = g, time = c(0, fit$time[ii]), estimate = c(0, fit$pstate[ii, j]),
               std.err = c(0, se), lower = c(0, lo), upper = c(0, hi),
               n.risk = c(if (length(ii)) fit$n.risk[ii][1L] else NA_real_, fit$n.risk[ii]),
               n.event = c(0, if (is.matrix(fit$n.event)) rowSums(fit$n.event[ii, , drop = FALSE]) else fit$n.event[ii]),
               n.censor = c(0, fit$n.censor[ii]), kind = "cif", stringsAsFactors = FALSE)
  })
  do.call(rbind, out)
}

.r4vn_surv_step_value <- function(curve, times) {
  vapply(times, function(tt) {
    ii <- which(curve$time <= tt)
    if (!length(ii)) return(curve$estimate[1L])
    curve$estimate[max(ii)]
  }, numeric(1))
}

.r4vn_surv_step_component <- function(curve, times, component, default = NA_real_) {
  vapply(times, function(tt) {
    ii <- which(curve$time <= tt)
    if (!length(ii)) return(default)
    curve[[component]][max(ii)]
  }, numeric(1))
}

.r4vn_surv_at_risk <- function(d, times, groups = NULL) {
  if (is.null(groups)) groups <- rep("All", nrow(d))
  lev <- unique(as.character(groups[!is.na(groups)]))
  out <- do.call(rbind, lapply(lev, function(g) {
    z <- d[as.character(groups) == g, , drop = FALSE]
    n <- vapply(times, function(tt) {
      if (".start" %in% names(z)) sum(z$.start <= tt & z$.time >= tt, na.rm = TRUE) else sum(z$.time >= tt, na.rm = TRUE)
    }, integer(1))
    data.frame(group = g, time = times, n.risk = n, stringsAsFactors = FALSE)
  }))
  rownames(out) <- NULL
  out
}

.r4vn_surv_km_at <- function(fit, times) {
  s <- summary(fit, times = times, extend = TRUE, data.frame = TRUE)
  if (!nrow(s)) return(data.frame())
  if (!"strata" %in% names(s)) s$strata <- "All"
  s$group <- .r4vn_surv_clean_strata(s$strata)
  s[, intersect(c("group", "time", "n.risk", "n.event", "n.censor", "surv", "std.err", "lower", "upper"), names(s)), drop = FALSE]
}

.r4vn_surv_cif_at <- function(curve, d, times) {
  groups <- unique(curve$group)
  risk <- .r4vn_surv_at_risk(d, times, if (".by" %in% names(d)) d$.by else NULL)
  out <- do.call(rbind, lapply(groups, function(g) {
    z <- curve[curve$group == g, , drop = FALSE]
    data.frame(group = g, time = times,
               cif = .r4vn_surv_step_value(z, times),
               std.err = .r4vn_surv_step_component(z, times, "std.err", NA_real_),
               lower = .r4vn_surv_step_component(z, times, "lower", 0),
               upper = .r4vn_surv_step_component(z, times, "upper", 0), stringsAsFactors = FALSE)
  }))
  merge(out, risk, by = c("group", "time"), all.x = TRUE, sort = FALSE)
}

.r4vn_surv_overview <- function(d, event_label = "Event") {
  n <- nrow(d)
  ev <- sum(d$.event == 1L, na.rm = TRUE)
  cens <- n - ev
  pt <- if (".start" %in% names(d)) sum(d$.time - d$.start, na.rm = TRUE) else sum(d$.time, na.rm = TRUE)
  data.frame(N = n, Events = ev, Censored = cens, `Event percent` = if (n) ev / n * 100 else NA_real_,
             `Person-time` = pt, check.names = FALSE)
}

.r4vn_surv_group_overview <- function(d) {
  if (!".by" %in% names(d)) return(NULL)
  groups <- levels(droplevels(d$.by))
  out <- lapply(groups, function(g) {
    z <- d[!is.na(d$.by) & d$.by == g, , drop = FALSE]
    ov <- .r4vn_surv_overview(z)
    data.frame(group = g, ov, stringsAsFactors = FALSE, check.names = FALSE)
  })
  do.call(rbind, out)
}

.r4vn_surv_auto_times <- function(time, n = 3L) {
  x <- as.numeric(time)
  x <- x[is.finite(x) & x > 0]
  if (!length(x)) return(NULL)
  probs <- seq(0, 1, length.out = n + 2L)[-c(1L, n + 2L)]
  ans <- as.numeric(stats::quantile(x, probs = probs, na.rm = TRUE, names = FALSE, type = 2))
  mx <- max(x)
  digits <- if (mx >= 20) 0L else if (mx >= 2) 1L else 2L
  ans <- sort(unique(round(ans, digits = digits)))
  ans[is.finite(ans) & ans > 0]
}

.r4vn_surv_interpretation <- function(x, alpha = 0.05) {
  rows <- list()
  add <- function(section, text) {
    if (!is.null(text) && length(text) && !is.na(text[1L]) && nzchar(text[1L])) {
      rows[[length(rows) + 1L]] <<- data.frame(
        section = section, interpretation = as.character(text[1L]),
        stringsAsFactors = FALSE
      )
    }
  }
  ov <- x$overview
  unit <- .r4vn_surv_unit_label(x$metadata$unit, TRUE)
  add("Overview", sprintf(
    "The analysis included %d observations, %d events (%.1f%%), and %d censored observations over %.2f person-%s.",
    ov$N[1L], ov$Events[1L], ov$`Event percent`[1L], ov$Censored[1L],
    ov$`Person-time`[1L], unit
  ))
  if (!is.null(x$followup) && nrow(x$followup) && is.finite(x$followup$median[1L])) {
    add("Follow-up", sprintf(
      "The reverse Kaplan-Meier median follow-up was %s %s.",
      .r4vn_surv_ci_text(x$followup$median[1L], x$followup$lower[1L],
                         x$followup$upper[1L], 2), unit
    ))
  }
  if (!is.null(x$logrank) && nrow(x$logrank) && is.finite(x$logrank$p[1L])) {
    add("Group comparison", sprintf(
      "The survival curves %s statistically different by the log-rank test (p %s).",
      if (x$logrank$p[1L] < alpha) "were" else "were not",
      .r4vn_surv_fmt_p(x$logrank$p[1L], 3)
    ))
  }
  if (!is.null(x$risk_compare) && nrow(x$risk_compare)) {
    z <- x$risk_compare[which.max(x$risk_compare$time), , drop = FALSE]
    parts <- character()
    if ("rr" %in% names(z) && is.finite(z$rr[1L])) {
      parts <- c(parts, paste0(
        "RR ", .r4vn_surv_ci_text(z$rr[1L], z$rr_lower[1L], z$rr_upper[1L], 2)
      ))
    }
    if ("rd" %in% names(z) && is.finite(z$rd[1L])) {
      parts <- c(parts, paste0(
        "RD ", .r4vn_surv_ci_text(z$rd[1L], z$rd_lower[1L], z$rd_upper[1L], 3)
      ))
    }
    if (length(parts)) add("Risk comparison", sprintf(
      "At time %.2f, %s versus %s: %s.", z$time[1L], z$comparison[1L],
      z$reference[1L], paste(parts, collapse = "; ")
    ))
  }
  if (!is.null(x$irr) && nrow(x$irr)) {
    z <- x$irr[if ("Overall" %in% x$irr$interval) which(x$irr$interval == "Overall")[1L] else 1L, , drop = FALSE]
    if (is.finite(z$irr[1L])) add("Incidence-rate comparison", sprintf(
      "%s versus %s: IRR %s, p %s.", z$comparison[1L], z$reference[1L],
      .r4vn_surv_ci_text(z$irr[1L], z$lower[1L], z$upper[1L], 2),
      .r4vn_surv_fmt_p(z$p[1L], 3)
    ))
  }
  model <- if (!is.null(x$cox$multi) && nrow(x$cox$multi)) x$cox$multi else x$cox$crude
  if (!is.null(model) && nrow(model)) {
    sig <- model[!model$reference & is.finite(model$p) & model$p < alpha, , drop = FALSE]
    if (nrow(sig)) {
      txt <- vapply(seq_len(nrow(sig)), function(i) sprintf(
        "%s (%s): HR %s, p %s",
        if ("variable_label" %in% names(sig)) sig$variable_label[i] else sig$variable[i],
        sig$level[i],
        .r4vn_surv_ci_text(sig$estimate[i], sig$lower[i], sig$upper[i], 2),
        .r4vn_surv_fmt_p(sig$p[i], 3)
      ), character(1))
      add("Cox model", paste0(
        "Statistically significant associations were observed for ",
        paste(txt, collapse = "; "),
        ". Hazard ratios describe associations and should not be interpreted as causal effects without an appropriate study design."
      ))
    } else {
      add("Cox model", "No modeled predictor had a two-sided p-value below 0.05 in the displayed Cox model.")
    }
  }
  if (!is.null(x$finegray$table) && nrow(x$finegray$table)) {
    fg <- x$finegray$table
    sig <- fg[!fg$reference & is.finite(fg$p) & fg$p < alpha, , drop = FALSE]
    if (nrow(sig)) {
      txt <- vapply(seq_len(nrow(sig)), function(i) sprintf(
        "%s (%s): SHR %s, p %s",
        if ("variable_label" %in% names(sig)) sig$variable_label[i] else sig$variable[i],
        sig$level[i],
        .r4vn_surv_ci_text(sig$estimate[i], sig$lower[i], sig$upper[i], 2),
        .r4vn_surv_fmt_p(sig$p[i], 3)
      ), character(1))
      add("Fine-Gray model", paste(txt, collapse = "; "))
    }
  }
  if (!is.null(x$ph) && nrow(x$ph) && "term" %in% names(x$ph) && "p" %in% names(x$ph)) {
    gi <- which(toupper(as.character(x$ph$term)) == "GLOBAL")
    if (length(gi) && is.finite(x$ph$p[gi[1L]])) {
      add("Model diagnostics", sprintf(
        "The global proportional-hazards test %s evidence against the PH assumption (p %s).",
        if (x$ph$p[gi[1L]] < alpha) "showed" else "did not show",
        .r4vn_surv_fmt_p(x$ph$p[gi[1L]], 3)
      ))
    }
  }
  if (isTRUE(x$metadata$competing)) {
    add("Competing risks", "Cumulative incidence was estimated with the Aalen-Johansen method; cause-specific Kaplan-Meier risk was not substituted for the competing-risk estimate.")
  }
  if (!is.null(x$rmst$difference) && nrow(x$rmst$difference)) {
    z <- x$rmst$difference[1L, , drop = FALSE]
    add("RMST comparison", sprintf(
      "%s minus %s: restricted mean survival-time difference %s, p %s.",
      z$comparison[1L], z$reference[1L],
      .r4vn_surv_ci_text(z$difference[1L], z$lower[1L], z$upper[1L], 2),
      .r4vn_surv_fmt_p(z$p[1L], 3)
    ))
  }
  if (!length(rows)) return(NULL)
  do.call(rbind, rows)
}

.r4vn_surv_reverse_followup <- function(d, ci = .95) {
  if (".start" %in% names(d)) return(NULL)
  fit <- survival::survfit(survival::Surv(.time, 1L - .event) ~ 1, data = d, conf.int = ci)
  tb <- summary(fit)$table
  if (is.matrix(tb)) tb <- tb[1L, ]
  getc <- function(pattern) {
    hit <- grep(pattern, names(tb), ignore.case = TRUE)
    if (length(hit)) unname(tb[hit[1L]]) else NA_real_
  }
  data.frame(median = getc("^median$"), lower = getc("LCL"), upper = getc("UCL"), stringsAsFactors = FALSE)
}

.r4vn_surv_median <- function(fit) {
  tb <- summary(fit)$table
  if (is.null(dim(tb))) tb <- matrix(tb, nrow = 1L, dimnames = list("All", names(tb)))
  cn <- colnames(tb)
  med <- grep("^median$", cn, ignore.case = TRUE)
  lo <- grep("LCL", cn, ignore.case = TRUE)
  hi <- grep("UCL", cn, ignore.case = TRUE)
  groups <- rownames(tb); if (is.null(groups)) groups <- "All"
  groups <- .r4vn_surv_clean_strata(groups)
  data.frame(group = groups,
             median = if (length(med)) tb[, med[1L]] else NA_real_,
             lower = if (length(lo)) tb[, lo[1L]] else NA_real_,
             upper = if (length(hi)) tb[, hi[1L]] else NA_real_, stringsAsFactors = FALSE)
}

.r4vn_surv_rmst <- function(fit, tau, ci = .95) {
  tb <- summary(fit, rmean = tau)$table
  if (is.null(dim(tb))) tb <- matrix(tb, nrow = 1L, dimnames = list("All", names(tb)))
  cn <- colnames(tb)
  im <- grep("rmean", cn, ignore.case = TRUE)
  ise <- grep("se\\(rmean\\)|se.*rmean", cn, ignore.case = TRUE)
  if (!length(im)) stop("RMST could not be extracted from `survfit`.", call. = FALSE)
  if (length(im) > 1L && length(ise)) im <- setdiff(im, ise)
  est <- tb[, im[1L]]
  se <- if (length(ise)) tb[, ise[1L]] else rep(NA_real_, length(est))
  z <- stats::qnorm(1 - (1 - ci) / 2)
  groups <- rownames(tb); if (is.null(groups)) groups <- "All"
  groups <- .r4vn_surv_clean_strata(groups)
  out <- data.frame(group = groups, tau = tau, rmst = est, se = se,
                    lower = est - z * se, upper = est + z * se, stringsAsFactors = FALSE)
  if (nrow(out) == 2L && all(is.finite(out$se))) {
    diff <- out$rmst[2L] - out$rmst[1L]
    sed <- sqrt(sum(out$se^2))
    out_diff <- data.frame(reference = out$group[1L], comparison = out$group[2L], difference = diff,
                           lower = diff - z * sed, upper = diff + z * sed,
                           p = 2 * stats::pnorm(-abs(diff / sed)), stringsAsFactors = FALSE)
  } else out_diff <- NULL
  list(table = out, difference = out_diff)
}

.r4vn_surv_person_time <- function(d, a = -Inf, b = Inf) {
  st <- if (".start" %in% names(d)) d$.start else rep(0, nrow(d))
  left <- pmax(st, a)
  right <- pmin(d$.time, b)
  sum(pmax(0, right - left), na.rm = TRUE)
}

.r4vn_surv_events_interval <- function(d, a = -Inf, b = Inf) {
  sum(d$.event == 1L & d$.time > a & d$.time <= b, na.rm = TRUE)
}

.r4vn_surv_rate_ci <- function(events, pt, scale = 100, ci = .95) {
  if (!is.finite(pt) || pt <= 0) return(c(rate = NA_real_, lower = NA_real_, upper = NA_real_))
  alpha <- 1 - ci
  rate <- events / pt * scale
  lo <- if (events == 0) 0 else 0.5 * stats::qchisq(alpha / 2, 2 * events) / pt * scale
  hi <- 0.5 * stats::qchisq(1 - alpha / 2, 2 * (events + 1)) / pt * scale
  c(rate = rate, lower = lo, upper = hi)
}

.r4vn_surv_rate_table <- function(d, at = NULL, mode = TRUE, scale = 100, ci = .95) {
  groups <- if (".by" %in% names(d)) as.character(d$.by) else rep("All", nrow(d))
  lev <- unique(groups[!is.na(groups)])
  mode2 <- if (is.character(mode)) tolower(mode[1L]) else "cumulative"
  out <- list(); k <- 0L
  for (g in lev) {
    z <- d[groups == g, , drop = FALSE]
    if (identical(mode2, "overall") || is.null(at) || !length(at)) {
      ranges <- list(c(-Inf, Inf)); labels <- "Overall"
    } else if (identical(mode2, "interval")) {
      edges <- c(0, sort(unique(at)))
      ranges <- lapply(seq_len(length(edges) - 1L), function(i) c(edges[i], edges[i + 1L]))
      labels <- paste0(edges[-length(edges)], "-", edges[-1L])
    } else if (identical(mode2, "all")) {
      aa <- sort(unique(at))
      edges <- c(0, aa)
      ranges <- c(
        list(c(-Inf, Inf)),
        lapply(aa, function(b) c(-Inf, b)),
        lapply(seq_len(length(edges) - 1L), function(i) c(edges[i], edges[i + 1L]))
      )
      labels <- c(
        "Overall", paste0("0-", aa),
        paste0(edges[-length(edges)], "-", edges[-1L], " (interval)")
      )
    } else {
      aa <- sort(unique(at)); ranges <- lapply(aa, function(b) c(-Inf, b)); labels <- paste0("0-", aa)
    }
    for (i in seq_along(ranges)) {
      a <- ranges[[i]][1L]; b <- ranges[[i]][2L]
      pt <- .r4vn_surv_person_time(z, a, b)
      ev <- .r4vn_surv_events_interval(z, a, b)
      rr <- .r4vn_surv_rate_ci(ev, pt, scale, ci)
      k <- k + 1L
      out[[k]] <- data.frame(group = g, interval = labels[i], start = if (is.finite(a)) a else NA_real_, end = if (is.finite(b)) b else NA_real_,
                             events = ev, person_time = pt, rate = rr["rate"], lower = rr["lower"], upper = rr["upper"], scale = scale,
                             stringsAsFactors = FALSE)
    }
  }
  do.call(rbind, out)
}

.r4vn_surv_compare_rates <- function(rate_table, ci = .95) {
  groups <- unique(rate_table$group)
  if (length(groups) != 2L) stop("`irr = TRUE` requires exactly two groups in `by`.", call. = FALSE)
  ints <- unique(rate_table$interval)
  out <- lapply(ints, function(intv) {
    z <- rate_table[rate_table$interval == intv, , drop = FALSE]
    z <- z[match(groups, z$group), , drop = FALSE]
    if (nrow(z) != 2L || any(z$person_time <= 0)) return(NULL)
    tst <- try(stats::poisson.test(z$events[c(2, 1)], T = z$person_time[c(2, 1)], conf.level = ci), silent = TRUE)
    if (inherits(tst, "try-error")) return(data.frame(interval = intv, reference = groups[1], comparison = groups[2], irr = NA_real_, lower = NA_real_, upper = NA_real_, p = NA_real_))
    est <- if (length(tst$estimate)) unname(tst$estimate[1L]) else (z$events[2] / z$person_time[2]) / (z$events[1] / z$person_time[1])
    data.frame(interval = intv, reference = groups[1], comparison = groups[2], irr = est,
               lower = tst$conf.int[1L], upper = tst$conf.int[2L], p = tst$p.value, stringsAsFactors = FALSE)
  })
  do.call(rbind, Filter(Negate(is.null), out))
}

.r4vn_surv_compare_risk <- function(risk_table, ci = .95, rr = FALSE, rd = FALSE) {
  groups <- unique(risk_table$group)
  if (length(groups) != 2L) stop("Risk comparisons require exactly two groups in `by`.", call. = FALSE)
  times <- unique(risk_table$time)
  zcrit <- stats::qnorm(1 - (1 - ci) / 2)
  out <- lapply(times, function(tt) {
    z <- risk_table[risk_table$time == tt, , drop = FALSE]
    z <- z[match(groups, z$group), , drop = FALSE]
    if (nrow(z) != 2L) return(NULL)
    r0 <- z$risk[1L]; r1 <- z$risk[2L]
    se0 <- z$se[1L]; se1 <- z$se[2L]
    ans <- data.frame(time = tt, reference = groups[1L], comparison = groups[2L], stringsAsFactors = FALSE)
    if (rr) {
      if (all(c(r0, r1, se0, se1) > 0, na.rm = FALSE)) {
        lrr <- log(r1 / r0); sel <- sqrt((se1 / r1)^2 + (se0 / r0)^2)
        ans$rr <- exp(lrr); ans$rr_lower <- exp(lrr - zcrit * sel); ans$rr_upper <- exp(lrr + zcrit * sel); ans$rr_p <- 2 * stats::pnorm(-abs(lrr / sel))
      } else ans[c("rr", "rr_lower", "rr_upper", "rr_p")] <- NA_real_
    }
    if (rd) {
      dif <- r1 - r0; sed <- sqrt(se1^2 + se0^2)
      ans$rd <- dif; ans$rd_lower <- dif - zcrit * sed; ans$rd_upper <- dif + zcrit * sed; ans$rd_p <- 2 * stats::pnorm(-abs(dif / sed))
    }
    ans
  })
  do.call(rbind, Filter(Negate(is.null), out))
}

.r4vn_surv_rhs <- function(map, extra = character(), interaction = NULL) {
  terms <- if (is.null(map) || !nrow(map)) character() else map$internal
  terms <- c(terms, extra)
  if (!is.null(interaction) && length(interaction) == 2L) terms <- c(terms, paste(interaction, collapse = "*"))
  if (!length(terms)) "1" else paste(unique(terms), collapse = " + ")
}

.r4vn_surv_fit_cox <- function(d, map, start = FALSE, strata_name = NULL, cluster_name = NULL,
                               frailty_name = NULL, interaction_internal = NULL, ties = "efron") {
  extra <- character()
  if (!is.null(strata_name)) extra <- c(extra, paste0("strata(", strata_name, ")"))
  if (!is.null(cluster_name)) extra <- c(extra, paste0("cluster(", cluster_name, ")"))
  if (!is.null(frailty_name)) extra <- c(extra, paste0("frailty(", frailty_name, ")"))
  rhs <- .r4vn_surv_rhs(map, extra, interaction_internal)
  f <- .r4vn_surv_formula(start, rhs)
  fit <- survival::coxph(f, data = d, ties = ties, x = TRUE, y = TRUE, model = TRUE)

  # Preserve the user-facing predictor names for postestimation. tabsurv() fits
  # Cox models on collision-safe internal columns (.x1, .x2, ...), but margins()
  # and predict() should accept the original variables (age, treatment, ...).
  if (!is.null(map) && nrow(map)) {
    pred <- d[, map$internal, drop = FALSE]
    names(pred) <- map$variable
    used_names <- rownames(fit$model)
    hit <- match(used_names, rownames(d))
    if (length(hit) && all(!is.na(hit))) pred <- pred[hit, , drop = FALSE]
    fit$.r4vn_prediction_data <- pred
    fit$.r4vn_prediction_map <- map[, intersect(c("variable", "internal", "type", "reference_index"), names(map)), drop = FALSE]
  }
  fit
}

.r4vn_surv_extract_cox <- function(fit, map, ci = .95) {
  if (is.null(map) || !nrow(map)) return(data.frame())
  sm <- summary(fit, conf.int = ci)
  cf <- sm$coefficients; cc <- sm$conf.int
  if (is.null(dim(cf))) cf <- matrix(cf, nrow = 1L, dimnames = list(names(stats::coef(fit))[1L], names(cf)))
  if (is.null(dim(cc))) cc <- matrix(cc, nrow = 1L, dimnames = list(rownames(cf)[1L], names(cc)))
  rn <- rownames(cf)
  pcol <- grep("Pr\\(", colnames(cf)); if (!length(pcol)) pcol <- ncol(cf)
  out <- list(); k <- 0L
  for (i in seq_len(nrow(map))) {
    sp <- map[i, , drop = FALSE]; internal <- sp$internal
    if (sp$type != "categorical") {
      hit <- which(rn == internal)
      k <- k + 1L
      if (!length(hit)) {
        out[[k]] <- data.frame(variable = sp$variable, level = "Per 1 unit", reference = FALSE,
                               estimate = NA_real_, lower = NA_real_, upper = NA_real_, p = NA_real_, stringsAsFactors = FALSE)
      } else {
        j <- hit[1L]
        out[[k]] <- data.frame(variable = sp$variable, level = "Per 1 unit", reference = FALSE,
                               estimate = cc[j, "exp(coef)"], lower = cc[j, grep("lower", colnames(cc), ignore.case = TRUE)[1L]],
                               upper = cc[j, grep("upper", colnames(cc), ignore.case = TRUE)[1L]], p = cf[j, pcol[1L]], stringsAsFactors = FALSE)
      }
    } else {
      f <- fit$model[[internal]]
      lev <- levels(f)
      coef_idx <- grep(paste0("^", internal), rn)
      nonref <- lev[-1L]
      for (j in seq_along(lev)) {
        k <- k + 1L
        if (j == 1L) {
          out[[k]] <- data.frame(variable = sp$variable, level = lev[j], reference = TRUE,
                                 estimate = 1, lower = 1, upper = 1, p = NA_real_, stringsAsFactors = FALSE)
        } else {
          pos <- j - 1L
          rowj <- if (pos <= length(coef_idx)) coef_idx[pos] else NA_integer_
          if (is.na(rowj)) {
            out[[k]] <- data.frame(variable = sp$variable, level = lev[j], reference = FALSE,
                                   estimate = NA_real_, lower = NA_real_, upper = NA_real_, p = NA_real_, stringsAsFactors = FALSE)
          } else {
            out[[k]] <- data.frame(variable = sp$variable, level = lev[j], reference = FALSE,
                                   estimate = cc[rowj, "exp(coef)"], lower = cc[rowj, grep("lower", colnames(cc), ignore.case = TRUE)[1L]],
                                   upper = cc[rowj, grep("upper", colnames(cc), ignore.case = TRUE)[1L]], p = cf[rowj, pcol[1L]], stringsAsFactors = FALSE)
          }
        }
      }
    }
  }
  ans <- do.call(rbind, out)
  labels <- stats::setNames(map$variable_label, map$variable)
  ans$variable_label <- unname(labels[ans$variable])
  ans[, c("variable", "variable_label", setdiff(names(ans), c("variable", "variable_label"))), drop = FALSE]
}


.r4vn_surv_extract_interaction <- function(fit, map, ci = .95) {
  if (is.null(fit)) return(NULL)
  sm <- summary(fit, conf.int = ci)
  cf <- sm$coefficients; cc <- sm$conf.int
  if (is.null(dim(cf))) cf <- matrix(cf, nrow = 1L, dimnames = list(names(stats::coef(fit))[1L], names(cf)))
  if (is.null(dim(cc))) cc <- matrix(cc, nrow = 1L, dimnames = list(rownames(cf)[1L], names(cc)))
  rn <- rownames(cf)
  hit <- grep(":", rn, fixed = TRUE)
  if (!length(hit)) return(NULL)
  pcol <- grep("Pr\\(", colnames(cf)); if (!length(pcol)) pcol <- ncol(cf)
  lo_col <- grep("lower", colnames(cc), ignore.case = TRUE)[1L]
  hi_col <- grep("upper", colnames(cc), ignore.case = TRUE)[1L]
  labels <- rn[hit]
  if (!is.null(map) && nrow(map)) {
    ord <- order(nchar(map$internal), decreasing = TRUE)
    for (i in ord) labels <- gsub(map$internal[i], map$variable[i], labels, fixed = TRUE)
  }
  data.frame(term = labels,
             estimate = cc[hit, "exp(coef)"], lower = cc[hit, lo_col], upper = cc[hit, hi_col],
             p = cf[hit, pcol[1L]], stringsAsFactors = FALSE)
}

.r4vn_surv_cox_diagnostics <- function(fit) {
  if (is.null(fit)) return(NULL)
  sm <- summary(fit)
  pulltest <- function(z, prefix) {
    if (is.null(z) || !length(z)) return(c(statistic = NA_real_, df = NA_real_, p = NA_real_))
    nm <- names(z)
    stat <- if ("test" %in% nm) z[["test"]] else z[1L]
    df <- if ("df" %in% nm) z[["df"]] else if (length(z) >= 2L) z[2L] else NA_real_
    pv <- if ("pvalue" %in% nm) z[["pvalue"]] else if (length(z) >= 3L) z[3L] else NA_real_
    c(statistic = unname(stat), df = unname(df), p = unname(pv))
  }
  lr <- pulltest(sm$logtest, "LR")
  wa <- pulltest(sm$waldtest, "Wald")
  sc <- pulltest(sm$sctest, "Score")
  conc <- sm$concordance
  data.frame(
    n = if (!is.null(fit$n)) fit$n else NA_real_,
    events = if (!is.null(fit$nevent)) fit$nevent else NA_real_,
    concordance = if (length(conc)) unname(conc[1L]) else NA_real_,
    concordance_se = if (length(conc) >= 2L) unname(conc[2L]) else NA_real_,
    AIC = tryCatch(stats::AIC(fit), error = function(e) NA_real_),
    LR_chisq = lr["statistic"], LR_df = lr["df"], LR_p = lr["p"],
    Wald_chisq = wa["statistic"], Wald_df = wa["df"], Wald_p = wa["p"],
    Score_chisq = sc["statistic"], Score_df = sc["df"], Score_p = sc["p"],
    stringsAsFactors = FALSE, check.names = FALSE
  )
}

.r4vn_surv_cox_models <- function(d, focal_map, all_map, adjusted_spec = NULL, multi_spec = NULL,
                                  start = FALSE, strata_name = NULL, cluster_name = NULL, frailty_name = NULL,
                                  ci = .95, ties = "efron", interaction_vars = NULL) {
  crude <- adjusted <- multi <- NULL
  crude_fits <- adjusted_fits <- list()
  interaction_table <- diagnostics <- NULL

  if (!is.null(focal_map) && nrow(focal_map)) {
    crude <- do.call(rbind, lapply(seq_len(nrow(focal_map)), function(i) {
      mp <- focal_map[i, , drop = FALSE]
      fit <- .r4vn_surv_fit_cox(d, mp, start, strata_name, cluster_name, frailty_name, ties = ties)
      crude_fits[[mp$variable]] <<- fit
      .r4vn_surv_extract_cox(fit, mp, ci)
    }))
  }

  if (!is.null(adjusted_spec)) {
    adjust_all <- isTRUE(adjusted_spec)
    if (adjust_all) adjvars <- focal_map$variable else adjvars <- adjusted_spec$variable
    adjusted <- do.call(rbind, lapply(seq_len(nrow(focal_map)), function(i) {
      focal <- focal_map[i, , drop = FALSE]
      covars <- setdiff(adjvars, focal$variable)
      use <- unique(c(focal$variable, covars))
      mp <- all_map[match(use, all_map$variable), , drop = FALSE]
      mp <- mp[!is.na(mp$variable), , drop = FALSE]
      fit <- .r4vn_surv_fit_cox(d, mp, start, strata_name, cluster_name, frailty_name, ties = ties)
      adjusted_fits[[focal$variable]] <<- fit
      .r4vn_surv_extract_cox(fit, focal, ci)
    }))
  }

  multi_fit <- NULL
  if (!is.null(multi_spec)) {
    if (isTRUE(multi_spec)) mp <- focal_map else mp <- all_map[all_map$variable %in% multi_spec$variable, , drop = FALSE]
    interaction_internal <- NULL
    if (!is.null(interaction_vars) && length(interaction_vars) == 2L) {
      interaction_internal <- mp$internal[match(interaction_vars, mp$variable)]
      if (anyNA(interaction_internal)) stop("Both interaction variables must be present in the multivariable model.", call. = FALSE)
    }
    multi_fit <- .r4vn_surv_fit_cox(d, mp, start, strata_name, cluster_name, frailty_name, interaction_internal, ties)
    multi <- .r4vn_surv_extract_cox(multi_fit, mp, ci)
    interaction_table <- .r4vn_surv_extract_interaction(multi_fit, mp, ci)
    diagnostics <- .r4vn_surv_cox_diagnostics(multi_fit)
  }

  list(crude = crude, adjusted = adjusted, multi = multi, interaction = interaction_table,
       diagnostics = diagnostics, crude_fits = crude_fits,
       adjusted_fits = adjusted_fits, multi_fit = multi_fit)
}

.r4vn_surv_ph <- function(fit) {
  if (is.null(fit)) return(NULL)
  z <- survival::cox.zph(fit, terms = TRUE, global = TRUE)
  tb <- as.data.frame(z$table)
  tb$term <- rownames(tb); rownames(tb) <- NULL
  names(tb)[names(tb) == "p"] <- "p"
  tb[, c("term", setdiff(names(tb), "term")), drop = FALSE]
}

.r4vn_surv_finegray <- function(d, map, start = FALSE, ci = .95, ties = "efron") {
  if (is.null(map) || !nrow(map)) stop("`finegray = TRUE` requires `vars` or `multi` covariates.", call. = FALSE)
  d$.fgid <- if (".id" %in% names(d)) d$.id else seq_len(nrow(d))
  lhs <- if (start) "survival::Surv(.start, .time, .status_ms)" else "survival::Surv(.time, .status_ms)"
  rhs <- paste(c(map$internal, ".fgid"), collapse = " + ")
  fgform <- stats::as.formula(paste(lhs, "~", rhs))
  fg <- survival::finegray(fgform, data = d, etype = "failure", id = if (start) d$.fgid else NULL)
  fitform <- stats::as.formula(paste("survival::Surv(fgstart, fgstop, fgstatus) ~", paste(map$internal, collapse = " + "), "+ cluster(.fgid)"))
  fit <- survival::coxph(fitform, data = fg, weights = fg[["fgwt"]], ties = ties, x = TRUE, y = TRUE, model = TRUE)
  list(fit = fit, table = .r4vn_surv_extract_cox(fit, map, ci), data = fg)
}

.r4vn_surv_ai_text <- function(x, digits = 2L, p_digits = 3L) {
  lines <- c(
    sprintf("Survival analysis: N=%d, events=%d, censored=%d.", x$overview$N[1L], x$overview$Events[1L], x$overview$Censored[1L])
  )
  if (!is.null(x$followup) && nrow(x$followup)) lines <- c(lines, sprintf("Median follow-up: %s.", .r4vn_surv_ci_text(x$followup$median[1], x$followup$lower[1], x$followup$upper[1], digits)))
  if (!is.null(x$logrank)) lines <- c(lines, sprintf("Log-rank p=%s.", .r4vn_surv_fmt_p(x$logrank$p[1], p_digits)))
  if (!is.null(x$cox$multi) && nrow(x$cox$multi)) {
    zz <- x$cox$multi[!x$cox$multi$reference, , drop = FALSE]
    lines <- c(lines, apply(zz, 1, function(r) sprintf("%s %s: HR %s, p=%s.", r[["variable"]], r[["level"]],
      .r4vn_surv_ci_text(as.numeric(r[["estimate"]]), as.numeric(r[["lower"]]), as.numeric(r[["upper"]]), digits),
      .r4vn_surv_fmt_p(as.numeric(r[["p"]]), p_digits))))
  }
  paste(lines, collapse = "\n")
}