tests/testthat/test-mixedbiastest.R

In mixedbiastest: Bias Diagnostic for Linear Mixed Models

## tests/testthat/test-mixedbiastest.R

## These tests are only run when NOT_CRAN = "true" because
## tests/testthat.R is wrapped in that condition.

test_that("mixedbiastest runs and returns expected components for simple model", {
  skip_if_not_installed("lme4")

  library(lme4)

  data("sleepstudy", package = "lme4")

  m1 <- lmer(Reaction ~ Days + (1 | Subject),
             data = sleepstudy,
             REML = FALSE)

  res <- mixedbiastest::mixedbiastest(m1, n_permutations = 40, verbose = FALSE)

  ## Basic structure
  expect_s3_class(res, "mixedbiastest")
  expect_true(is.list(res))
  expect_true(all(c("results_table", "permutation_values", "model") %in% names(res)))

  ## Results table structure
  rt <- res$results_table
  expect_true(is.data.frame(rt))
  expect_true(all(c("Fixed_Effect", "Mean_Permuted_Value",
                    "P_Value", "Bias_Estimate") %in% names(rt)))
  expect_gt(nrow(rt), 0L)

  ## permutation_values list
  pv <- res$permutation_values
  expect_true(is.list(pv))
  expect_equal(names(pv), rt$Fixed_Effect)
})

test_that("mixedbiastest works with multiple diagonal random-intercept terms", {
  skip_if_not_installed("lme4")

  library(lme4)

  set.seed(123)
  n_grp1 <- 10
  n_grp2 <- 8
  n_obs  <- 200

  grp1   <- factor(sample(seq_len(n_grp1), size = n_obs, replace = TRUE))
  grp2   <- factor(sample(seq_len(n_grp2), size = n_obs, replace = TRUE))
  x      <- runif(n_obs)

  u1 <- rnorm(n_grp1, 0, 1)
  u2 <- rnorm(n_grp2, 0, 0.5)
  e  <- rnorm(n_obs, 0, 0.3)
  y  <- 1 + 2 * x + u1[grp1] + u2[grp2] + e

  dat <- data.frame(y = y, x = x, grp1 = grp1, grp2 = grp2)

  m <- lmer(y ~ x + (1 | grp1) + (1 | grp2),
            data = dat, REML = FALSE)

  res <- mixedbiastest::mixedbiastest(m, n_permutations = 30, verbose = FALSE)

  expect_s3_class(res, "mixedbiastest")
  ## Should have as many rows as fixed effects
  expect_equal(nrow(res$results_table), length(lme4::fixef(m)))
})

test_that("mixedbiastest accepts separate diagonal random intercept and slope terms", {
  skip_if_not_installed("lme4")

  library(lme4)

  data("sleepstudy", package = "lme4")

  ## (1|Subject) + (0+Days|Subject) -> separate 1x1 VarCorr blocks -> diagonal G
  m_sep <- lmer(Reaction ~ Days + (1 | Subject) + (0 + Days | Subject),
                data = sleepstudy, REML = FALSE)

  res <- mixedbiastest::mixedbiastest(m_sep, n_permutations = 20, verbose = FALSE)

  expect_s3_class(res, "mixedbiastest")
  expect_equal(nrow(res$results_table), length(lme4::fixef(m_sep)))
})

test_that("mixedbiastest rejects models with non-diagonal G (correlated intercept & slope)", {
  skip_if_not_installed("lme4")

  library(lme4)

  data("sleepstudy", package = "lme4")

  ## (Days | Subject) -> 2x2 VarCorr with nonzero off-diagonals
  m_nd <- lmer(Reaction ~ Days + (Days | Subject),
               data = sleepstudy, REML = FALSE)

  expect_error(
    mixedbiastest::mixedbiastest(m_nd, n_permutations = 10, verbose = FALSE),
    "not diagonal",
    ignore.case = TRUE
  )
})

test_that("mixedbiastest with non-unit weights is currently allowed (no special handling)", {
  skip_if_not_installed("lme4")

  library(lme4)

  data("sleepstudy", package = "lme4")

  set.seed(2025)
  w <- runif(nrow(sleepstudy), min = 0.5, max = 1.5)

  m_w <- lmer(Reaction ~ Days + (1 | Subject),
              data = sleepstudy,
              weights = w,
              REML = FALSE)

  ## NOTE: At present, mixedbiastest does not explicitly check or forbid weights.
  ## This test just ensures it *runs*; the documentation should make clear
  ## that the method assumes R = sigma^2 I, so weighted fits are at user's risk.
  res <- mixedbiastest::mixedbiastest(m_w, n_permutations = 10, verbose = FALSE)
  expect_s3_class(res, "mixedbiastest")
})

test_that("plot.mixedbiastest returns a ggplot object", {
  skip_if_not_installed("lme4")
  skip_if_not_installed("ggplot2")

  library(lme4)
  library(ggplot2)

  data("sleepstudy", package = "lme4")

  m1 <- lmer(Reaction ~ Days + (1 | Subject),
             data = sleepstudy, REML = FALSE)

  res <- mixedbiastest::mixedbiastest(m1, n_permutations = 30, verbose = FALSE)

  p <- plot(res)

  expect_s3_class(p, "ggplot")
})

test_that("list_fixed_effects returns correct names and indices", {
  skip_if_not_installed("lme4")

  library(lme4)

  data("sleepstudy", package = "lme4")

  m1 <- lmer(Reaction ~ Days + (1 | Subject),
             data = sleepstudy, REML = FALSE)

  lf <- mixedbiastest::list_fixed_effects(m1)

  expect_true(is.data.frame(lf))
  expect_true(all(c("Index", "Coefficient") %in% names(lf)))
  expect_equal(lf$Index, seq_len(nrow(lf)))
  expect_equal(lf$Coefficient, names(lme4::fixef(m1)))
})

test_that("print.mixedbiastest prints without error and returns object invisibly", {
  skip_if_not_installed("lme4")

  library(lme4)

  data("sleepstudy", package = "lme4")

  m1 <- lmer(Reaction ~ Days + (1 | Subject),
             data = sleepstudy, REML = FALSE)

  res <- mixedbiastest::mixedbiastest(m1, n_permutations = 20, verbose = FALSE)

  ## We expect printed output, so we *don't* use expect_silent().
  ## Instead we check that the return value is invisible and equal to res.
  expect_invisible(out <- print(res))
  expect_identical(out, res)
})


test_that("mixedbiastest works on the example_data dataset", {
  skip_if_not_installed("lme4")

  library(lme4)

  ## example_data is shipped with the mixedbiastest package
  data("example_data", package = "mixedbiastest")

  m_ex <- lmer(y ~ x + (1 | group),
               data = example_data,
               REML = FALSE)

  res <- mixedbiastest::mixedbiastest(m_ex, n_permutations = 40, verbose = FALSE)

  expect_s3_class(res, "mixedbiastest")
  expect_equal(nrow(res$results_table), length(lme4::fixef(m_ex)))
})