tests/testthat/test-r4vn-1-4-1-upgrade.R

In R4VN: Health Data Analysis and Publication-Ready Reporting

# Regression tests for the R4VN 1.4.1 usability/statistics upgrade.

test_that("genvar supports quick vectors, compact repetition, and row expansion", {
  usedf(clear = TRUE, quiet = TRUE)
  on.exit(usedf(clear = TRUE, quiet = TRUE), add = TRUE)

  genvar(weight = c(29, 26, 13, 23, 23))
  d <- usedf(quiet = TRUE)
  expect_equal(nrow(d), 5L)
  expect_equal(d$weight, c(29, 26, 13, 23, 23))

  genvar(smoking = c(1, 0), times = c(4, 4))
  d <- usedf(quiet = TRUE)
  expect_equal(nrow(d), 8L)
  expect_equal(d$smoking, c(rep(1, 4), rep(0, 4)))
  expect_equal(d$weight[1:5], c(29, 26, 13, 23, 23))
  expect_true(all(is.na(d$weight[6:8])))

  genvar(group = c("A", "B"), each = 5)
  d <- usedf(quiet = TRUE)
  expect_equal(nrow(d), 10L)
  expect_equal(d$group, rep(c("A", "B"), each = 5))
  expect_true(all(is.na(d$weight[9:10])))
})

test_that("normality tests accept many variables and hierarchical by", {
  set.seed(1401)
  d <- data.frame(
    x = rnorm(80), y = rnorm(80, 2),
    province = factor(rep(c("A", "B"), each = 40)),
    sex = factor(rep(rep(c("F", "M"), each = 20), 2))
  )
  z <- normtest(vars = vars(x, y), by = vars(province, sex), data = d,
                method = c("shapiro", "jarque.bera"), show = FALSE)
  expect_s3_class(z, "r4vn_stat")
  expect_true(all(c("Variable", "Stratum", "Group", "Method") %in% names(z$sections$Test)))
  expect_equal(nrow(z$sections$Test), 16L)

  s <- swilk(vars = vars(x, y), by = vars(province, sex), data = d, show = FALSE)
  expect_s3_class(s, "r4vn_stat")
  expect_true(all(s$sections$Test$Method == "Shapiro-Wilk"))
})

test_that("t test, ANOVA, and Kruskal-Wallis expose requested extensions", {
  set.seed(1402)
  d <- data.frame(
    y = c(rnorm(30, 0), rnorm(30, .5), rnorm(30, 1.2)),
    group = factor(rep(c("A", "B", "C"), each = 30)),
    province = factor(rep(c("North", "South"), length.out = 90))
  )
  d$group2 <- factor(ifelse(d$group == "A", "A", "Other"))

  tt <- ttest(y, by = group2, data = d, effect = TRUE, show = FALSE)
  expect_true("Effect size" %in% names(tt$sections))

  av <- anova(y, by = group, data = d, posthoc = "tukey", show = FALSE)
  expect_true("Effect size" %in% names(av$sections))
  expect_true("Tukey honestly significant difference test" %in% names(av$sections))

  avb <- anova(y, by = group, data = d, posthoc = "bonferroni", show = FALSE)
  expect_true("Bonferroni-adjusted pairwise t-tests" %in% names(avb$sections))
  avb_posthoc <- avb$sections[["Bonferroni-adjusted pairwise t-tests"]]
  expect_true("p.adjusted" %in% names(avb_posthoc))
  expect_identical(attr(avb_posthoc, "adjust"), "bonferroni")

  avi <- anovai(c(25, 10, 2), c(25, 12, 2.5), c(25, 15, 3),
                group.names = c("A", "B", "C"), posthoc = "games-howell", show = FALSE)
  expect_true("Games-Howell multiple comparisons" %in% names(avi$sections))

  avib <- anovai(c(25, 10, 2), c(25, 12, 2.5), c(25, 15, 3),
                 group.names = c("A", "B", "C"), posthoc = "bonferroni", show = FALSE)
  expect_true("Bonferroni-adjusted pairwise t-tests" %in% names(avib$sections))

  kw <- kwallis(y, by = group, data = d, posthoc = "dunn", effect = TRUE, show = FALSE)
  expect_true("Dunn's test (Holm-adjusted p-values)" %in% names(kw$sections))
  expect_true("Effect size" %in% names(kw$sections))

  h <- anova(y, by = vars(province, group), data = d, posthoc = "pairwise", show = FALSE)
  expect_s3_class(h, "r4vn_stat")
  expect_true("Stratum" %in% names(h$sections[[1L]]))
})

test_that("hierarchical by convention works in proportion tests and tab", {
  set.seed(1403)
  d <- data.frame(
    event = rbinom(120, 1, .35),
    sex = factor(rep(c("F", "M"), 60)),
    province = factor(rep(c("A", "B", "C"), each = 40))
  )
  attr(d$province, "label") <- "Province"
  attr(d$sex, "label") <- "Sex"
  p <- prtest(event, by = vars(province, sex), data = d, event = 1, show = FALSE)
  expect_s3_class(p, "r4vn_stat")
  expect_true(any(vapply(p$sections, function(x) is.data.frame(x) && "Stratum" %in% names(x), logical(1))))

  tb <- tab(vars = vars(event), by = vars(province, sex), data = d, show = FALSE)
  expect_s3_class(tb, "r4vn_tab")
  expect_equal(tb$hierarchical_by$by, "sex")
  expect_equal(tb$hierarchical_by$strata, "province")
  expect_equal(tb$hierarchical_by$by_label, "Sex")
  expect_true(all(grepl("^Province \\(.*\\)$", tb$superby_levels)))
})

test_that("graph commands accept vars, hierarchical by, and reference lines", {
  d <- data.frame(
    x = 1:24,
    y1 = seq(10, 33), y2 = seq(20, 43),
    province = factor(rep(c("A", "B"), each = 12)),
    sex = factor(rep(c("F", "M"), 12))
  )
  attr(d$province, "label") <- "Province"
  attr(d$sex, "label") <- "Sex"
  attr(d$y1, "label") <- "Outcome 1"
  attr(d$y2, "label") <- "Outcome 2"
  f <- tempfile(fileext = ".png")
  h <- ghist(data = d, vars = vars(y1, y2), by = vars(province, sex),
             normal = TRUE, xline = 25, ref_lty = 3, combine = TRUE,
             ncol = 2, file = f, show = FALSE)
  expect_s3_class(h, "r4vn_graph_set")
  expect_equal(length(h$graphs), 8L)
  expect_true(h$combined)
  expect_equal(h$ncol, 2L)
  expect_true(file.exists(f) && file.info(f)$size > 0)
  expect_true(all(grepl("Outcome [12] \\| Province \\([AB]\\) > Sex \\([FM]\\)", h$labels)))
  expect_silent(print(h))
  expect_silent(print(h$graphs[[1L]]))

  s <- gscatter(data = d, x = x, vars = vars(y1, y2), by = vars(province, sex),
                yline = 25, xline = 12, show = FALSE)
  expect_s3_class(s, "r4vn_graph_set")
  expect_equal(length(s$graphs), 4L)

  lf <- tempfile(fileext = ".png")
  l <- gline(data = d, x = x, vars = vars(y1, y2),
             by = vars(province, sex), combine = TRUE, ncol = 2,
             file = lf, show = FALSE)
  expect_s3_class(l, "r4vn_graph_set")
  expect_true(file.exists(lf) && file.info(lf)$size > 0)

  # Deprecated names remain compatibility aliases, but cannot be mixed with
  # their canonical replacements.
  expect_s3_class(ghist(d, y1, vline = 25, show = FALSE), "r4vn_graph")
  expect_error(ghist(d, y1, xline = 25, vline = 25, show = FALSE), "Use only `xline`")
})

test_that("varform returns all nine Tukey ladder transformations", {
  d <- data.frame(x = c(1, 2, 3, 5, 8, 13, 21, 34))
  z <- varform(x, data = d, show = FALSE)
  expect_s3_class(z, "r4vn_stat")
  expect_equal(nrow(z$sections[[1L]]), 9L)
  expect_equal(length(z$raw$transformations[[1L]]), 9L)
})

test_that("active-model margins, predict, and lincom interoperate", {
  set.seed(1404)
  d <- data.frame(
    y = rbinom(160, 1, .4),
    age = rnorm(160, 45, 10),
    htn = rbinom(160, 1, .3)
  )
  usedf(d, quiet = TRUE)
  on.exit(usedf(clear = TRUE, quiet = TRUE), add = TRUE)

  m <- logistic(y, c.age, i.htn, data = d, event = 1, show = FALSE)
  mg <- margins(at = at(age = c(35, 45, 55), htn = c(0, 1)), show = FALSE)
  expect_s3_class(mg, "r4vn_margins")
  expect_equal(nrow(mg$raw$margins), 6L)
  expect_true(all(mg$raw$margins$Margin >= 0 & mg$raw$margins$Margin <= 1))

  out <- predict(newvar = phat, type = "probability", show = FALSE)
  expect_true("phat" %in% names(out))
  expect_equal(length(out$phat), nrow(d))

  rr <- predict(m$raw$model, newdata = d[1:5, ], type = "response")
  expect_equal(length(rr), 5L)

  cn <- names(stats::coef(m$raw$model))
  nonint <- setdiff(cn, "(Intercept)")
  expect_true(length(nonint) >= 1L)
  lc <- lincom(paste0("`", nonint[1L], "` + 1"), model = m, show = FALSE)
  expect_s3_class(lc, "r4vn_stat")
})

test_that("binary Poisson exposes event and risk-ratio output", {
  set.seed(1405)
  d <- data.frame(y = rbinom(180, 1, .3), x = rnorm(180), g = factor(rep(c("A", "B"), 90)))
  z <- poisson(y, c.x, i.g, data = d, event = 1, rr = TRUE, vce = "robust", show = FALSE)
  expect_s3_class(z, "r4vn_stat")
  expect_true(isTRUE(z$raw$binary))
  expect_equal(z$raw$event, "1")
  expect_true("Risk ratios" %in% names(z$sections))
})

test_that("nptrend supports numeric and binary ordered trends", {
  d <- data.frame(
    dose = ordered(rep(c("Low", "Mid", "High"), each = 30), levels = c("Low", "Mid", "High")),
    score = c(rnorm(30, 10), rnorm(30, 12), rnorm(30, 14)),
    event = c(rbinom(30, 1, .1), rbinom(30, 1, .3), rbinom(30, 1, .6))
  )
  a <- nptrend(score, by = dose, data = d, method = "cuzick", show = FALSE)
  b <- nptrend(event, by = dose, data = d, method = "cochran-armitage", event = 1, show = FALSE)
  expect_s3_class(a, "r4vn_stat")
  expect_s3_class(b, "r4vn_stat")
})

test_that("nlregress fits spline shapes and becomes postestimation-ready", {
  set.seed(1406)
  d <- data.frame(x = seq(0, 10, length.out = 100))
  d$y <- 2 + .5 * d$x - .08 * d$x^2 + rnorm(100, 0, .5)
  z <- nlregress(y, x, data = d, spline = "natural", df = 4, show = FALSE)
  expect_s3_class(z, "r4vn_stat")
  expect_s3_class(z$raw$model, "lm")
  expect_equal(z$raw$x, "x")
})

test_that("qregress works when optional quantreg is installed", {
  skip_if_not_installed("quantreg")
  set.seed(1407)
  d <- data.frame(y = rnorm(80), x = rnorm(80))
  z <- qregress(y, c.x, data = d, tau = c(.25, .5, .75), show = FALSE)
  expect_s3_class(z, "r4vn_stat")
  expect_equal(length(z$raw$models), 3L)
})


test_that("epi reports crude-versus-MH OR differences and hierarchical by", {
  # Construct two provinces, each with two sex strata and nondegenerate 2x2 data.
  d <- expand.grid(province = c("A", "B"), sex = c("F", "M"),
                   exposure = c(0, 1), outcome = c(0, 1), rep = 1:8,
                   KEEP.OUT.ATTRS = FALSE, stringsAsFactors = FALSE)
  # Tilt the event distribution so estimates are finite and not all identical.
  keep <- !(d$exposure == 0 & d$outcome == 1 & d$rep > 3) &
          !(d$exposure == 1 & d$outcome == 0 & d$rep > 5)
  d <- d[keep, ]
  d$province <- factor(d$province); d$sex <- factor(d$sex)

  one <- epi(outcome, exposure, by = sex, data = d[d$province == "A", ],
             event = 1, exposed = 1, show = FALSE)
  adj <- one$raw$stratified
  expect_true(all(c("or.percent.diff.mh", "or.percent.diff.crude") %in% names(adj)))

  h <- epi(outcome, exposure, by = vars(province, sex), data = d,
           event = 1, exposed = 1, show = FALSE)
  expect_s3_class(h, "r4vn_stat")
  expect_equal(h$hierarchical_by$by, "sex")
  expect_equal(h$hierarchical_by$strata, "province")
})

test_that("haven application-control errors are identified for fallback guidance", {
  msg <- R4VN:::.r4vn_haven_block_message("LoadLibrary failure: An Application Control policy has blocked this file")
  expect_match(msg, "Application Control")
  expect_match(msg, "readstata13")
})

test_that("saved vars selectors retain hierarchical by semantics", {
  d <- data.frame(
    y = rnorm(72),
    province = factor(rep(c("A", "B"), each = 36)),
    sex = factor(rep(rep(c("F", "M"), each = 18), 2))
  )
  g <- vars(province, sex)
  z <- normtest(y, by = g, data = d, method = "shapiro", show = FALSE)
  expect_s3_class(z, "r4vn_stat")
  expect_true(all(c("Stratum", "Group") %in% names(z$sections$Test)))
  expect_equal(sort(unique(z$sections$Test$Stratum)), c("province (A)", "province (B)"))

  # y is continuous; R4VN therefore uses the explicit c. declaration.
  tb <- tab(vars = vars(c.y), by = g, data = d, show = FALSE)
  expect_equal(tb$hierarchical_by$by, "sex")
  expect_equal(tb$hierarchical_by$strata, "province")

  gr <- ghist(data = d, x = y, by = g, show = FALSE)
  expect_s3_class(gr, "r4vn_graph_set")
  expect_equal(length(gr$graphs), 4L)
})



test_that("tab bounds Fisher computation for large sparse categorical tables", {
  d <- data.frame(
    idcat = seq_len(30),
    group = factor(rep(c("A", "B"), 15))
  )
  # Deliberately leave idcat unprefixed to exercise the categorical safety
  # guard. This is not the recommended analysis; c.idcat is recommended for a
  # genuinely continuous variable. The important contract is that tab()
  # returns instead of entering an unbounded exact Fisher calculation.
  z <- tab(vars = vars(idcat), by = group, data = d, show = FALSE)
  expect_s3_class(z, "r4vn_tab")
})

test_that("tabsurv accepts hierarchical by selectors", {
  skip_if_not_installed("survival")
  set.seed(1408)
  d <- data.frame(
    time = rexp(120, .1),
    event = rbinom(120, 1, .55),
    province = factor(rep(c("A", "B"), each = 60)),
    sex = factor(rep(c("F", "M"), 60))
  )
  z <- tabsurv(time, event, by = vars(province, sex), data = d, show = FALSE)
  expect_s3_class(z, "r4vn_surv")
  expect_equal(z$hierarchical_by$by, "sex")
  expect_equal(z$hierarchical_by$strata, "province")
  expect_equal(length(z$hierarchical_results), 2L)
  expect_true(all(vapply(z$hierarchical_results, inherits, logical(1), what = "r4vn_surv")))
})

test_that("spline postestimation varies the original predictor", {
  set.seed(1409)
  d <- data.frame(x = seq(1, 10, length.out = 120))
  d$y <- 1 + .7 * d$x - .05 * d$x^2 + rnorm(120, 0, .4)
  nlregress(y, x, data = d, spline = "natural", df = 4, show = FALSE)
  mg <- margins(at = at(x = c(2, 5, 8)), show = FALSE)
  expect_s3_class(mg, "r4vn_margins")
  expect_equal(mg$raw$margins$x, c(2, 5, 8))
  expect_true(all(is.finite(mg$raw$margins$Margin)))
})

test_that("Cox models become active for margins and prediction", {
  skip_if_not_installed("survival")
  set.seed(1410)
  d <- data.frame(
    time = rexp(160, .08), event = rbinom(160, 1, .65),
    age = rnorm(160, 50, 10)
  )
  usedf(d, quiet = TRUE)
  on.exit(usedf(clear = TRUE, quiet = TRUE), add = TRUE)
  cox(time, event, vars = vars(c.age), data = d, failure = 1, ph = FALSE, show = FALSE)
  mg <- margins(at = at(age = c(40, 50, 60)), type = "response", show = FALSE)
  expect_s3_class(mg, "r4vn_margins")
  expect_true(all(mg$raw$margins$Margin > 0))
  pp <- predict(newvar = rrisk, type = "risk", show = FALSE)
  expect_true("rrisk" %in% names(pp))
  expect_equal(length(pp$rrisk), nrow(d))
})