tests/testthat/test-models-lrtest.R

In R4VN: Health Data Analysis and Publication-Ready Reporting

# Tests for compact regression syntax and likelihood-ratio model comparison

test_that("compact logistic syntax supports interaction terms", {
  set.seed(20260820)

  n <- 500
  d <- data.frame(
    age = rnorm(n, 45, 12),
    occupation = factor(
      sample(c("Office", "Worker", "Other"), n, TRUE),
      levels = c("Office", "Worker", "Other")
    ),
    intervention = factor(
      sample(c("No", "Yes"), n, TRUE),
      levels = c("No", "Yes")
    )
  )

  eta <- -1.5 +
    0.02 * d$age +
    0.30 * (d$occupation == "Worker") +
    0.20 * (d$occupation == "Other") +
    0.40 * (d$intervention == "Yes") +
    0.60 * (d$occupation == "Worker" & d$intervention == "Yes")

  d$outcome <- factor(
    rbinom(n, 1, plogis(eta)),
    levels = 0:1,
    labels = c("No", "Yes")
  )

  m1 <- logistic(
    outcome,
    c.age,
    i.occupation,
    i.intervention,
    data = d,
    event = "Yes",
    show = FALSE
  )

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

  expect_s3_class(m1, "r4vn_stat")
  expect_s3_class(m2, "r4vn_stat")

  tl <- attr(stats::terms(m2$raw$model), "term.labels")
  expect_true("occupation:intervention" %in% tl)

  z <- lrtest(m1, m2, show = FALSE)
  expect_s3_class(z, "r4vn_stat")
  expect_equal(nrow(z$raw$table), 2L)
  expect_equal(z$raw$table$df[2L], 2L)
  expect_true(is.finite(z$raw$table$LR[2L]))
  expect_true(z$raw$table$p[2L] >= 0 && z$raw$table$p[2L] <= 1)
})


test_that("inline vars syntax is captured without changing public vars", {
  set.seed(20260821)

  n <- 300
  d <- data.frame(
    age = rnorm(n, 40, 10),
    occupation = factor(
      sample(c("Office", "Worker", "Other"), n, TRUE),
      levels = c("Office", "Worker", "Other")
    ),
    intervention = factor(
      sample(c("No", "Yes"), n, TRUE),
      levels = c("No", "Yes")
    )
  )

  d$outcome <- factor(
    rbinom(n, 1, 0.35),
    levels = 0:1,
    labels = c("No", "Yes")
  )

  m <- logistic(
    outcome,
    vars = vars(
      c.age,
      ib2.occupation,
      i.intervention,
      ib2.occupation*i.intervention
    ),
    data = d,
    event = "Yes",
    show = FALSE
  )

  expect_s3_class(m, "r4vn_stat")
  mf <- stats::model.frame(m$raw$model)
  expect_true(is.factor(mf$occupation))
  expect_identical(levels(mf$occupation)[1L], "Worker")

  tl <- attr(stats::terms(m$raw$model), "term.labels")
  expect_true("occupation:intervention" %in% tl)
})


test_that("ib prefix prefers literal level and otherwise uses positive index fallback", {
  set.seed(20260822)

  n <- 240

  d1 <- data.frame(
    group = factor(sample(c("1", "2", "3"), n, TRUE), levels = c("1", "2", "3")),
    x = rnorm(n)
  )
  d1$y <- factor(
    rbinom(n, 1, 0.4),
    levels = 0:1,
    labels = c("No", "Yes")
  )

  m1 <- logistic(
    y,
    ib2.group,
    c.x,
    data = d1,
    event = "Yes",
    show = FALSE
  )

  expect_identical(
    levels(stats::model.frame(m1$raw$model)$group)[1L],
    "2"
  )

  d2 <- d1
  d2$group <- factor(
    c("Low", "Middle", "High")[as.integer(d1$group)],
    levels = c("Low", "Middle", "High")
  )

  m2 <- logistic(
    y,
    ib2.group,
    c.x,
    data = d2,
    event = "Yes",
    show = FALSE
  )

  expect_identical(
    levels(stats::model.frame(m2$raw$model)$group)[1L],
    "Middle"
  )
})


test_that("b prefix keeps existing R4VN factor-position semantics", {
  set.seed(20260823)

  n <- 220
  d <- data.frame(
    group = factor(
      sample(c("A", "B", "C"), n, TRUE),
      levels = c("A", "B", "C")
    ),
    x = rnorm(n)
  )

  d$y <- factor(
    rbinom(n, 1, 0.35),
    levels = 0:1,
    labels = c("No", "Yes")
  )

  m <- logistic(
    y,
    b3.group,
    c.x,
    data = d,
    event = "Yes",
    show = FALSE
  )

  expect_identical(
    levels(stats::model.frame(m$raw$model)$group)[1L],
    "C"
  )
})


test_that("three logistic models are compared sequentially", {
  set.seed(20260824)

  n <- 400
  d <- data.frame(
    age = rnorm(n, 45, 12),
    sex = factor(sample(c("Female", "Male"), n, TRUE)),
    treatment = factor(sample(c("No", "Yes"), n, TRUE))
  )

  pr <- plogis(
    -1.5 +
      0.025 * d$age +
      0.35 * (d$sex == "Male") +
      0.45 * (d$treatment == "Yes")
  )

  d$y <- factor(
    rbinom(n, 1, pr),
    levels = 0:1,
    labels = c("No", "Yes")
  )

  m1 <- logistic(
    y,
    c.age,
    data = d,
    event = "Yes",
    show = FALSE
  )

  m2 <- logistic(
    y,
    c.age,
    i.sex,
    data = d,
    event = "Yes",
    show = FALSE
  )

  m3 <- logistic(
    y,
    c.age,
    i.sex,
    i.treatment,
    data = d,
    event = "Yes",
    show = FALSE
  )

  z <- lrtest(m1, m2, m3, show = FALSE)

  expect_equal(nrow(z$raw$table), 3L)
  expect_true(is.na(z$raw$table$p[1L]))
  expect_true(is.finite(z$raw$table$p[2L]))
  expect_true(is.finite(z$raw$table$p[3L]))
})


test_that("Poisson compact syntax and lrtest work with the same exposure", {
  set.seed(20260825)

  n <- 450
  d <- data.frame(
    time = runif(n, 0.5, 5),
    age = rnorm(n, 45, 12),
    sex = factor(sample(c("Female", "Male"), n, TRUE)),
    treatment = factor(sample(c("No", "Yes"), n, TRUE))
  )

  rate <- exp(
    -0.4 +
      0.01 * d$age +
      0.25 * (d$sex == "Male") +
      0.35 * (d$treatment == "Yes")
  )

  d$cases <- rpois(n, lambda = d$time * rate)

  p1 <- poisson(
    cases,
    c.age,
    i.sex,
    i.treatment,
    data = d,
    exposure = time,
    show = FALSE
  )

  p2 <- poisson(
    cases,
    c.age,
    i.sex*i.treatment,
    data = d,
    exposure = time,
    show = FALSE
  )

  z <- lrtest(p1, p2, show = FALSE)

  expect_s3_class(z, "r4vn_stat")
  expect_identical(z$raw$type, "Poisson regression")
  expect_true(is.finite(z$raw$table$LR[2L]))
})


test_that("lrtest rejects different analytic samples", {
  set.seed(20260826)

  n <- 250
  d <- data.frame(
    age = rnorm(n),
    bmi = rnorm(n)
  )
  d$bmi[1:25] <- NA_real_
  d$y <- factor(
    rbinom(n, 1, 0.4),
    levels = 0:1,
    labels = c("No", "Yes")
  )

  m1 <- logistic(
    y,
    c.age,
    data = d,
    event = "Yes",
    show = FALSE
  )

  m2 <- logistic(
    y,
    c.age,
    c.bmi,
    data = d,
    event = "Yes",
    show = FALSE
  )

  expect_error(
    lrtest(m1, m2, show = FALSE),
    "different analytic samples|different observations"
  )
})


test_that("lrtest rejects robust R4VN fits", {
  set.seed(20260827)

  n <- 240
  d <- data.frame(
    x = rnorm(n),
    z = rnorm(n)
  )
  d$y <- factor(
    rbinom(n, 1, 0.4),
    levels = 0:1,
    labels = c("No", "Yes")
  )

  m1 <- logistic(
    y,
    c.x,
    data = d,
    event = "Yes",
    vce = "robust",
    show = FALSE
  )

  m2 <- logistic(
    y,
    c.x,
    c.z,
    data = d,
    event = "Yes",
    vce = "robust",
    show = FALSE
  )

  expect_error(
    lrtest(m1, m2, show = FALSE),
    "vce"
  )
})


test_that("native glm objects are accepted and quasi models are rejected", {
  set.seed(20260828)

  n <- 250
  d <- data.frame(
    x = rnorm(n),
    z = rnorm(n)
  )

  d$y <- rbinom(
    n,
    1,
    plogis(-1 + 0.5 * d$x + 0.4 * d$z)
  )

  g1 <- stats::glm(
    y ~ x,
    data = d,
    family = stats::binomial()
  )

  g2 <- stats::glm(
    y ~ x + z,
    data = d,
    family = stats::binomial()
  )

  expect_s3_class(
    lrtest(g1, g2, show = FALSE),
    "r4vn_stat"
  )

  d$count <- rpois(n, exp(0.2 + 0.2 * d$x))

  q1 <- stats::glm(
    count ~ x,
    data = d,
    family = stats::quasipoisson()
  )

  q2 <- stats::glm(
    count ~ x + z,
    data = d,
    family = stats::quasipoisson()
  )

  expect_error(
    lrtest(q1, q2, show = FALSE),
    "quasi"
  )
})


test_that("ordinary linear regression is directed to the nested F test", {
  set.seed(20260829)

  d <- data.frame(
    y = rnorm(200),
    x = rnorm(200),
    z = rnorm(200)
  )

  m1 <- stats::lm(y ~ x, data = d)
  m2 <- stats::lm(y ~ x + z, data = d)

  expect_error(
    lrtest(m1, m2, show = FALSE),
    "linear regression"
  )
})


test_that("negative binomial models are accepted when MASS is available", {
  skip_if_not_installed("MASS")

  set.seed(20260830)

  n <- 300
  d <- data.frame(
    x = rnorm(n),
    z = rnorm(n)
  )

  mu <- exp(0.3 + 0.25 * d$x + 0.2 * d$z)
  d$y <- MASS::rnegbin(n, mu = mu, theta = 2)

  m1 <- MASS::glm.nb(y ~ x, data = d)
  m2 <- MASS::glm.nb(y ~ x + z, data = d)

  z <- lrtest(m1, m2, show = FALSE)
  expect_identical(z$raw$type, "Negative binomial regression")
})


test_that("Cox models are accepted when survival is available", {
  skip_if_not_installed("survival")

  set.seed(20260831)

  n <- 300
  d <- data.frame(
    x = rnorm(n),
    z = rnorm(n)
  )

  rate <- exp(0.25 * d$x + 0.20 * d$z)
  event_time <- rexp(n, rate = rate)
  censor_time <- rexp(n, rate = 0.5)

  d$time <- pmin(event_time, censor_time)
  d$status <- as.integer(event_time <= censor_time)

  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
  )

  z <- lrtest(m1, m2, show = FALSE)
  expect_identical(z$raw$type, "Cox regression")
})