tests/testthat/test-tabsurv.R

In R4VN: Health Data Analysis and Publication-Ready Reporting

test_that("tabsurv returns KM, risk and incidence rates", {
  skip_if_not_installed("survival")
  set.seed(20260808)
  n <- 240
  d <- data.frame(
    trt = factor(sample(c("Control", "Treatment"), n, TRUE)),
    age = round(rnorm(n, 55, 12)),
    sex = factor(sample(c("Female", "Male"), n, TRUE))
  )
  haz <- 0.08 * exp(-0.45 * (d$trt == "Treatment") + 0.02 * (d$age - 55))
  te <- rexp(n, haz); tc <- runif(n, 1, 8)
  d$time <- pmin(te, tc)
  d$death <- as.integer(te <= tc)

  z <- tabsurv(time, death, by = trt, data = d, failure = 1,
               at = c(1, 3, 5), risk = TRUE, rate = TRUE,
               show = FALSE)
  expect_s3_class(z, "r4vn_surv")
  expect_equal(z$overview$N, n)
  expect_true(nrow(z$risk) == 6L)
  expect_true(all(c("rate", "lower", "upper") %in% names(z$rate)))
  expect_true(is.data.frame(z$logrank))
})

test_that("missing Cox covariates do not reduce KM sample", {
  skip_if_not_installed("survival")
  set.seed(20260808)
  n <- 150
  d <- data.frame(time = rexp(n, .15), event = rbinom(n, 1, .6),
                  age = rnorm(n, 50, 10), sex = factor(sample(c("F", "M"), n, TRUE)))
  d$age[1:25] <- NA
  z <- tabsurv(time, event, data = d, vars = vars(c.age, sex),
               cox = TRUE, multi = TRUE, show = FALSE)
  expect_equal(z$overview$N, n)
  expect_true(z$cox$diagnostics$n <= n - 25)
})

test_that("b2 reference is respected in Cox model", {
  skip_if_not_installed("survival")
  set.seed(20260808)
  n <- 180
  d <- data.frame(time = rexp(n, .1), event = rbinom(n, 1, .65),
                  sex = factor(sample(c("Female", "Male"), n, TRUE), levels = c("Female", "Male")))
  z <- tabsurv(time, event, data = d, vars = vars(b2.sex), cox = TRUE, show = FALSE)
  expect_equal(z$cox$crude$level[z$cox$crude$reference], "Male")
})

test_that("competing risks produce Aalen-Johansen CIF", {
  skip_if_not_installed("survival")
  set.seed(20260808)
  n <- 220
  t1 <- rexp(n, .08); t2 <- rexp(n, .05); tc <- runif(n, 2, 10)
  tm <- pmin(t1, t2, tc)
  status <- ifelse(tm == t1, 1L, ifelse(tm == t2, 2L, 0L))
  d <- data.frame(time = tm, status = status,
                  trt = factor(sample(c("A", "B"), n, TRUE)), age = rnorm(n, 60, 8))
  z <- tabsurv(time, status, data = d, failure = 1, compete = 2,
               by = trt, at = c(2, 5), risk = TRUE, show = FALSE)
  expect_true(isTRUE(z$metadata$competing))
  expect_true(nrow(z$risk) == 4L)
  expect_true(all(z$risk$risk >= 0 & z$risk$risk <= 1, na.rm = TRUE))
})

test_that("Fine-Gray and PH outputs are reusable", {
  skip_if_not_installed("survival")
  set.seed(20260808)
  n <- 260
  trt <- factor(sample(c("A", "B"), n, TRUE))
  age <- rnorm(n, 58, 10)
  t1 <- rexp(n, .07 * exp(.02 * (age - 58)))
  t2 <- rexp(n, .04); tc <- runif(n, 2, 12)
  tm <- pmin(t1, t2, tc)
  status <- ifelse(tm == t1, 1L, ifelse(tm == t2, 2L, 0L))
  d <- data.frame(time = tm, status = status, trt = trt, age = age)
  z <- tabsurv(time, status, data = d, failure = 1, compete = 2,
               vars = vars(c.age, trt), finegray = TRUE, show = FALSE)
  expect_true(is.data.frame(z$finegray$table))
})

test_that("Andersen-Gill recurrent event syntax accepts start-stop-id", {
  skip_if_not_installed("survival")
  d <- data.frame(
    id = rep(1:40, each = 2),
    start = rep(c(0, 1), 40),
    stop = rep(c(1, 2), 40),
    event = rbinom(80, 1, .25),
    trt = factor(rep(sample(c("A", "B"), 40, TRUE), each = 2)),
    age = rep(rnorm(40, 55, 8), each = 2)
  )
  z <- tabsurv(stop, event, data = d, start = start, id = id,
               vars = vars(trt, c.age), recurrent = "ag", show = FALSE)
  expect_s3_class(z, "r4vn_surv")
  expect_true(!is.null(z$cox$multi_fit))
})

test_that("automatic report provides a complete two-group survival report", {
  skip_if_not_installed("survival")
  set.seed(20260901)
  n <- 260
  d <- data.frame(
    trt = factor(rep(c("Control", "Treatment"), each = n / 2)),
    age = rnorm(n, 56, 9)
  )
  event_time <- rexp(n, .10 * exp(-.40 * (d$trt == "Treatment")))
  censor_time <- runif(n, 2, 12)
  d$time <- pmin(event_time, censor_time)
  d$event <- as.integer(event_time <= censor_time)

  z <- tabsurv(time, event, by = trt, data = d, show = FALSE)

  expect_identical(z$metadata$report, "auto")
  expect_true(is.data.frame(z$descriptive))
  expect_equal(nrow(z$descriptive), 2L)
  expect_true(is.data.frame(z$risk))
  expect_true(is.data.frame(z$rate))
  expect_true(all(z$rate$interval == "Overall"))
  expect_true(is.data.frame(z$logrank))
  expect_true(is.data.frame(z$risk_compare))
  expect_true(is.data.frame(z$irr))
  expect_true(is.list(z$rmst))
  expect_null(z$interpretation)
  expect_null(z$lifetable)
  expect_true(is.list(z$tables))
  expect_true(is.list(z$plots))
})

test_that("interpretation and life table are explicit opt-in modules", {
  skip_if_not_installed("survival")
  set.seed(20260906)
  n <- 180
  d <- data.frame(
    time = rexp(n, .11),
    event = rbinom(n, 1, .62),
    trt = factor(rep(c("A", "B"), each = n / 2))
  )

  z0 <- tabsurv(time, event, by = trt, data = d,
                report = "custom", show = FALSE)
  expect_null(z0$interpretation)
  expect_null(z0$lifetable)
  expect_identical(formals(tabsurv)$interpretation, FALSE)
  expect_identical(formals(tabsurv)$lifetable, FALSE)

  z <- tabsurv(time, event, by = trt, data = d,
               report = "custom", lifetable = TRUE,
               interpretation = TRUE, show = FALSE)
  expect_true(is.data.frame(z$lifetable))
  expect_true(nrow(z$lifetable) > 0L)
  expect_true(all(c(
    "group", "time", "n.risk", "n.event", "n.censor",
    "conditional_survival", "survival", "cumulative_risk",
    "std.error", "lower", "upper"
  ) %in% names(z$lifetable)))
  expect_equal(sum(z$lifetable$n.event), sum(d$event))
  expect_identical(z$estimates$life_table, z$lifetable)
  expect_true("Life table" %in% names(z$tables))
  expect_true(is.data.frame(z$interpretation))
})

test_that("cuminc reports cumulative incidence at exactly requested times", {
  skip_if_not_installed("survival")
  set.seed(20260907)
  n <- 220
  d <- data.frame(
    time = rexp(n, .08),
    event = rbinom(n, 1, .67),
    trt = factor(rep(c("A", "B"), each = n / 2))
  )

  z <- tabsurv(time, event, by = trt, data = d,
               cuminc = c(6, 12, 24), report = "custom", show = FALSE)
  expect_true(is.data.frame(z$cuminc))
  expect_identical(z$cuminc, z$risk)
  expect_equal(sort(unique(z$cuminc$time)), c(6, 12, 24))
  expect_equal(nrow(z$cuminc), 6L)
  expect_true(all(z$cuminc$risk >= 0 & z$cuminc$risk <= 1, na.rm = TRUE))
  expect_identical(z$metadata$cuminc, c(6, 12, 24))
})

test_that("cuminc uses Aalen-Johansen CIF with competing risks", {
  skip_if_not_installed("survival")
  set.seed(20260908)
  n <- 260
  t1 <- rexp(n, .07)
  t2 <- rexp(n, .05)
  tc <- runif(n, 4, 30)
  tm <- pmin(t1, t2, tc)
  status <- ifelse(tm == t1, 1L, ifelse(tm == t2, 2L, 0L))
  d <- data.frame(time = tm, status = status,
                  trt = factor(rep(c("A", "B"), each = n / 2)))

  z <- tabsurv(time, status, by = trt, data = d,
               failure = 1, compete = 2,
               cuminc = c(6, 12, 24), report = "custom", show = FALSE)
  expect_true(isTRUE(z$metadata$competing))
  expect_equal(sort(unique(z$cuminc$time)), c(6, 12, 24))
  expect_true(all(z$cuminc$risk >= 0 & z$cuminc$risk <= 1, na.rm = TRUE))
})

test_that("custom report preserves the concise R4VN 1.5 defaults", {
  skip_if_not_installed("survival")
  set.seed(20260902)
  d <- data.frame(
    time = rexp(120, .12), event = rbinom(120, 1, .6),
    trt = factor(rep(c("A", "B"), 60))
  )
  z <- tabsurv(time, event, by = trt, data = d,
               report = "custom", show = FALSE)

  expect_null(z$risk)
  expect_null(z$rate)
  expect_null(z$rmst)
  expect_null(z$graph)
  expect_true(is.data.frame(z$logrank))
})

test_that("automatic report activates Cox models and diagnostics with vars", {
  skip_if_not_installed("survival")
  set.seed(20260903)
  n <- 300
  d <- data.frame(
    age = rnorm(n, 58, 8),
    sex = factor(sample(c("Female", "Male"), n, TRUE)),
    trt = factor(sample(c("Control", "Treatment"), n, TRUE))
  )
  event_time <- rexp(n, .08 * exp(.025 * (d$age - 58) - .35 * (d$trt == "Treatment")))
  censor_time <- runif(n, 3, 14)
  d$time <- pmin(event_time, censor_time)
  d$event <- as.integer(event_time <= censor_time)
  attr(d$age, "label") <- "Age (years)"
  attr(d$trt, "label") <- "Treatment group"

  z <- tabsurv(time, event, by = trt,
               vars = vars(c.age, sex, trt), data = d, show = FALSE)

  expect_true(is.data.frame(z$cox$crude))
  expect_true("variable_label" %in% names(z$cox$crude))
  expect_true("Age (years)" %in% z$cox$crude$variable_label)
  expect_true(is.data.frame(z$cox$multi))
  expect_true(is.data.frame(z$cox$diagnostics))
  expect_true(is.data.frame(z$ph))
  expect_true(is.list(z$models))
  expect_true(is.list(z$diagnostics))
})

test_that("full report includes overall, cumulative and interval incidence rates", {
  skip_if_not_installed("survival")
  set.seed(20260904)
  d <- data.frame(
    time = rexp(180, .15), event = rbinom(180, 1, .65),
    trt = factor(rep(c("A", "B"), 90))
  )
  z <- tabsurv(time, event, by = trt, data = d,
               report = "full", show = FALSE)

  expect_true(any(z$rate$interval == "Overall"))
  expect_true(any(grepl("interval", z$rate$interval, fixed = TRUE)))
  expect_true(any(grepl("^0-", z$rate$interval)))
})

test_that("hierarchical by uses the same complete report contract", {
  skip_if_not_installed("survival")
  set.seed(20260905)
  n <- 320
  d <- data.frame(
    site = factor(rep(c("North", "South"), each = n / 2)),
    trt = factor(rep(rep(c("A", "B"), each = n / 4), 2)),
    time = rexp(n, .10), event = rbinom(n, 1, .6)
  )
  z <- tabsurv(time, event, by = vars(site, trt), data = d, show = FALSE)

  expect_s3_class(z, "r4vn_surv_hierarchical")
  expect_equal(length(z$hierarchical_results), 2L)
  expect_true(is.data.frame(z$descriptive))
  expect_true(is.list(z$tables))
  expect_true(is.list(z$plots))
  expect_false(exists(".r4vn_tabsurv_legacy", envir = asNamespace("R4VN"), inherits = FALSE))
})