tests/testthat/test-lrtest-model-syntax.R

In R4VN: Health Data Analysis and Publication-Ready Reporting

testthat::test_that("compact logistic syntax supports factor declarations and interactions", {
  set.seed(20260820)
  n <- 500L
  d <- data.frame(
    age = stats::rnorm(n, 45, 12),
    occupation = factor(sample(c("Office", "Clinical", "Other"), n, TRUE)),
    treatment = factor(sample(c("Control", "Intervention"), n, TRUE))
  )
  lp <- -2 + 0.025 * d$age +
    0.35 * (d$occupation == "Clinical") +
    0.25 * (d$occupation == "Other") +
    0.45 * (d$treatment == "Intervention") +
    0.55 * (d$occupation == "Clinical") * (d$treatment == "Intervention")
  d$outcome <- factor(stats::rbinom(n, 1, stats::plogis(lp)), levels = 0:1, labels = c("No", "Yes"))

  m1 <- R4VN::logistic(outcome, c.age, i.occupation, i.treatment,
                       data = d, event = "Yes", show = FALSE)
  m2 <- R4VN::logistic(outcome, c.age, i.occupation*i.treatment,
                       data = d, event = "Yes", show = FALSE)

  testthat::expect_s3_class(m1, "r4vn_stat")
  testthat::expect_s3_class(m2, "r4vn_stat")
  testthat::expect_true(any(grepl(":", attr(stats::terms(m2$raw$model), "term.labels"), fixed = TRUE)))

  z <- R4VN::lrtest(m1, m2, show = FALSE)
  testthat::expect_s3_class(z, "r4vn_stat")
  testthat::expect_equal(z$raw$comparison$LR.df[2L], 2L)
  testthat::expect_true(is.finite(z$raw$comparison$LR.chi2[2L]))
})

testthat::test_that("inline vars syntax supports ib reference and interaction", {
  set.seed(20260821)
  n <- 350L
  d <- data.frame(
    age = stats::rnorm(n, 50, 10),
    occupation = factor(sample(c("A", "B", "C"), n, TRUE), levels = c("A", "B", "C")),
    treatment = factor(sample(c("No", "Yes"), n, TRUE), levels = c("No", "Yes"))
  )
  d$outcome <- factor(stats::rbinom(n, 1, 0.35), levels = 0:1, labels = c("No", "Yes"))

  fit <- R4VN::logistic(
    outcome,
    vars = R4VN::vars(c.age, ib2.occupation*i.treatment),
    data = d,
    event = "Yes",
    show = FALSE
  )

  testthat::expect_s3_class(fit, "r4vn_stat")
  mf <- stats::model.frame(fit$raw$model)
  ref_col <- grep("r4vn_factor_index", names(mf), value = TRUE)
  testthat::expect_true(length(ref_col) >= 1L)
  testthat::expect_identical(levels(mf[[ref_col[1L]]])[1L], "B")
})

testthat::test_that("legacy b2 syntax remains usable in regression models", {
  set.seed(20260822)
  n <- 250L
  d <- data.frame(
    group = factor(sample(c("A", "B", "C"), n, TRUE), levels = c("A", "B", "C")),
    y = factor(stats::rbinom(n, 1, .4), levels = 0:1, labels = c("No", "Yes"))
  )
  fit <- R4VN::logistic(y, b2.group, data = d, event = "Yes", show = FALSE)
  mf <- stats::model.frame(fit$raw$model)
  ref_col <- grep("r4vn_factor_index", names(mf), value = TRUE)
  testthat::expect_identical(levels(mf[[ref_col[1L]]])[1L], "B")
})

testthat::test_that("lrtest supports sequential three-model comparisons", {
  set.seed(20260823)
  n <- 450L
  d <- data.frame(
    age = stats::rnorm(n),
    sex = factor(sample(c("F", "M"), n, TRUE)),
    trt = factor(sample(c("No", "Yes"), n, TRUE))
  )
  d$y <- factor(stats::rbinom(n, 1, stats::plogis(-.4 + .2*d$age)), levels = 0:1, labels = c("No", "Yes"))

  m1 <- R4VN::logistic(y, c.age, data = d, event = "Yes", show = FALSE)
  m2 <- R4VN::logistic(y, c.age, i.sex, data = d, event = "Yes", show = FALSE)
  m3 <- R4VN::logistic(y, c.age, i.sex, i.trt, data = d, event = "Yes", show = FALSE)
  z <- R4VN::lrtest(m1, m2, m3, show = FALSE)

  testthat::expect_equal(nrow(z$raw$comparison), 3L)
  testthat::expect_true(all(is.finite(z$raw$comparison$LR.chi2[-1L])))
  testthat::expect_true(all(z$raw$comparison$LR.df[-1L] > 0L))
})

testthat::test_that("lrtest works with Poisson regression", {
  set.seed(20260824)
  n <- 500L
  d <- data.frame(
    age = stats::rnorm(n),
    group = factor(sample(c("A", "B", "C"), n, TRUE)),
    time = stats::runif(n, .5, 3)
  )
  mu <- exp(.3 + .15*d$age + .25*(d$group == "B")) * d$time
  d$cases <- stats::rpois(n, mu)

  m1 <- R4VN::poisson(cases, c.age, data = d, exposure = time, show = FALSE)
  m2 <- R4VN::poisson(cases, c.age, i.group, data = d, exposure = time, show = FALSE)
  z <- R4VN::lrtest(m1, m2, show = FALSE)

  testthat::expect_s3_class(z, "r4vn_stat")
  testthat::expect_equal(z$raw$comparison$LR.df[2L], 2L)
})

testthat::test_that("lrtest accepts native GLM objects", {
  set.seed(20260825)
  n <- 300L
  d <- data.frame(x = stats::rnorm(n), z = stats::rnorm(n))
  d$y <- stats::rbinom(n, 1, stats::plogis(-.5 + .4*d$x))
  m1 <- stats::glm(y ~ x, family = stats::binomial(), data = d)
  m2 <- stats::glm(y ~ x + z, family = stats::binomial(), data = d)
  out <- R4VN::lrtest(m1, m2, show = FALSE)
  testthat::expect_s3_class(out, "r4vn_stat")
  testthat::expect_equal(out$raw$comparison$LR.df[2L], 1L)
})

testthat::test_that("lrtest rejects models fitted to different samples", {
  set.seed(20260826)
  d <- data.frame(y = stats::rbinom(100, 1, .4), x = stats::rnorm(100), z = stats::rnorm(100))
  d$z[1:10] <- NA_real_
  m1 <- stats::glm(y ~ x, family = stats::binomial(), data = d)
  m2 <- stats::glm(y ~ x + z, family = stats::binomial(), data = d)
  testthat::expect_error(R4VN::lrtest(m1, m2, show = FALSE), "different numbers of observations|different observations")
})

testthat::test_that("lrtest rejects quasi likelihood", {
  set.seed(20260827)
  d <- data.frame(y = stats::rpois(150, 2), x = stats::rnorm(150), z = stats::rnorm(150))
  m1 <- stats::glm(y ~ x, family = stats::quasipoisson(), data = d)
  m2 <- stats::glm(y ~ x + z, family = stats::quasipoisson(), data = d)
  testthat::expect_error(R4VN::lrtest(m1, m2, show = FALSE), "quasi")
})

testthat::test_that("lrtest rejects ordinary linear regression", {
  d <- data.frame(y = 1:30, x = stats::rnorm(30), z = stats::rnorm(30))
  m1 <- stats::lm(y ~ x, data = d)
  m2 <- stats::lm(y ~ x + z, data = d)
  testthat::expect_error(R4VN::lrtest(m1, m2, show = FALSE), "F test")
})

testthat::test_that("negative binomial models are supported when MASS is installed", {
  testthat::skip_if_not_installed("MASS")
  set.seed(20260828)
  n <- 300L
  d <- data.frame(x = stats::rnorm(n), z = stats::rnorm(n))
  d$y <- MASS::rnegbin(n, mu = exp(.5 + .3*d$x), theta = 1.5)
  m1 <- MASS::glm.nb(y ~ x, data = d)
  m2 <- MASS::glm.nb(y ~ x + z, data = d)
  out <- R4VN::lrtest(m1, m2, show = FALSE)
  testthat::expect_s3_class(out, "r4vn_stat")
})

testthat::test_that("Cox models are supported when survival is installed", {
  testthat::skip_if_not_installed("survival")
  set.seed(20260829)
  n <- 250L
  d <- data.frame(time = stats::rexp(n), status = stats::rbinom(n, 1, .7), x = stats::rnorm(n), z = stats::rnorm(n))
  m1 <- survival::coxph(survival::Surv(time, status) ~ x, data = d, model = TRUE, x = TRUE)
  m2 <- survival::coxph(survival::Surv(time, status) ~ x + z, data = d, model = TRUE, x = TRUE)
  out <- R4VN::lrtest(m1, m2, show = FALSE)
  testthat::expect_s3_class(out, "r4vn_stat")
})

testthat::test_that("redundant plain main effects are absorbed by starred categorical interaction", {
  set.seed(20260830)
  n <- 320L
  d <- data.frame(
    occupation = factor(sample(c("A", "B", "C"), n, TRUE)),
    preterm = factor(sample(c("No", "Yes"), n, TRUE))
  )
  d$outcome <- factor(stats::rbinom(n, 1, .35), levels = 0:1, labels = c("No", "Yes"))

  fit <- R4VN::logistic(
    outcome, occupation, preterm, i.occupation*i.preterm,
    data = d, event = "Yes", show = FALSE
  )

  mm <- stats::model.matrix(fit$raw$model)
  testthat::expect_false(any(is.na(stats::coef(fit$raw$model))))
  testthat::expect_equal(qr(mm)$rank, ncol(mm))
})