tests/testthat/test-glmmTMB.R

In modelbased: Estimation of Model-Based Predictions, Contrasts and Means

skip_on_cran()
skip_if_not_installed("glmmTMB")
skip_if_not_installed("emmeans")
skip_if_not_installed("marginaleffects", minimum_version = "0.29.0")

data <- glmmTMB::Salamanders
model <- suppressWarnings(glmmTMB::glmmTMB(
  count ~ mined + (1 | site),
  ziformula = ~mined,
  family = poisson,
  data = data
))

model2 <- suppressWarnings(glmmTMB::glmmTMB(
  count ~ cover + mined + (1 | site),
  ziformula = ~ cover + mined,
  family = glmmTMB::truncated_poisson,
  data = data
))


test_that("estimate_means - glmmTMB", {
  ## emmeans for zero-inflated model, count component
  estim <- suppressMessages(estimate_means(model, backend = "emmeans"))
  estim2 <- as.data.frame(emmeans::emmeans(model, ~mined, type = "response"))
  expect_equal(estim$Rate, estim2$rate, tolerance = 1e-3)

  ## emmeans for zero-inflated model, zero-inflated component
  estim <- suppressMessages(estimate_means(model, component = "zi", backend = "emmeans"))
  estim2 <- as.data.frame(emmeans::emmeans(model, ~mined, component = "zi", type = "response"))
  expect_equal(estim$Rate, estim2$rate, tolerance = 1e-3)

  ## marginaleffects for zero-inflated model, count component
  estim1 <- estimate_means(model, by = "mined", backend = "marginaleffects", predict = "inverse_link")
  estim2 <- suppressMessages(estimate_means(model, backend = "marginaleffects", predict = "inverse_link"))
  expect_equal(estim1$Mean, estim2$Mean, tolerance = 1e-3)
  estim3 <- estimate_means(model, by = "mined", backend = "emmeans")
  expect_equal(estim1$Mean, estim3$Rate, tolerance = 1e-1)

  ## marginaleffects for zero-inflated model, zero-inflation probabilities
  estim1 <- estimate_means(model, by = "mined", backend = "marginaleffects", predict = "zprob")
  estim2 <- estimate_means(model, backend = "marginaleffects", predict = "zprob")
  expect_equal(estim1$Probability, estim2$Probability, tolerance = 1e-3)
  # validated against ggeffects
  expect_equal(estim1$Probability, c(0.75755, 0.35508), tolerance = 1e-3)
  estim_me <- suppressWarnings(marginaleffects::avg_predictions(
    model,
    newdata = insight::get_datagrid(model, by = "mined", factors = "all", include_random = TRUE),
    by = "mined",
    re.form = NULL,
    type = "zprob"
  ))
  expect_equal(estim2$Probability, estim_me$estimate, tolerance = 1e-3)
})


test_that("estimate_contrasts - glmmTMB", {
  ## contrasts emmeans for zero-inflated model, count component
  estim1 <- suppressMessages(estimate_contrasts(model, backend = "emmeans"))
  expect_identical(dim(estim1), c(1L, 9L))
  # validated against ggeffects
  expect_equal(estim1$Difference, -2.32874, tolerance = 1e-3)
  expect_identical(c(estim1$Level1[1], estim1$Level2[1]), c("yes", "no"))

  ## contrasts marginaleffects for zero-inflated model, count component
  estim3 <- suppressMessages(estimate_contrasts(model, contrast = "mined", backend = "marginaleffects"))
  expect_equal(estim3$Difference, 1.99344, tolerance = 1e-3)
  estim_me <- suppressWarnings(marginaleffects::avg_predictions(
    model,
    newdata = insight::get_datagrid(model, by = "mined", factors = "all", include_random = TRUE),
    by = "mined",
    hypothesis = ~pairwise,
    re.form = NULL
  ))
  expect_equal(estim3$Difference, estim_me$estimate, tolerance = 1e-3)

  # select default for contrast automatically works
  estim4 <- suppressMessages(estimate_contrasts(model, backend = "marginaleffects"))
  expect_equal(estim3$Difference, estim4$Difference, tolerance = 1e-3)

  ## contrasts emmeans for zero-inflated model, zero-inflation probability component
  estim1 <- suppressMessages(estimate_contrasts(model, component = "zi", backend = "emmeans"))
  estim2 <- predict(model, type = "zprob", newdata = insight::get_datagrid(model, "mined"))
  expect_identical(dim(estim1), c(1L, 9L))
  expect_equal(estim1$Difference * -1, diff(estim2), tolerance = 1e-1)
  expect_identical(c(estim1$Level1[1], estim1$Level2[1]), c("yes", "no"))

  ## contrasts marginaleffects for zero-inflated model, zero-inflation probability component
  estim3 <- suppressMessages(estimate_contrasts(model, predict = "zprob", backend = "marginaleffects"))
  expect_equal(estim3$Difference, diff(estim2), tolerance = 1e-1)
})


test_that("estimate_slope - glmmTMB", {
  skip_on_os("mac")
  estim <- suppressMessages(estimate_slopes(model2, trend = "cover", by = "mined", regrid = "response", backend = "emmeans")) # nolint
  estim2 <- as.data.frame(emmeans::emtrends(model2, "mined", var = "cover", regrid = "response"))
  expect_equal(estim$Slope, estim2$cover.trend, tolerance = 1e-2)
  estim1 <- estimate_slopes(model2, trend = "cover", by = "mined", backend = "marginaleffects")
  datagrid <- insight::get_datagrid(model2, by = "mined")
  estim2 <- suppressWarnings(marginaleffects::avg_slopes(model2, newdata = datagrid, variables = "cover", by = "mined"))
  expect_equal(estim1$Slope, estim2$estimate, tolerance = 1e-2)
})


test_that("estimate_expectation - glmmTMB", {
  estim <- suppressMessages(estimate_expectation(model2))
  expect_identical(dim(estim), c(nrow(data), 8L))
})


test_that("estimate_link - glmmTMB", {
  estim <- suppressMessages(estimate_link(model2, preserve_range = FALSE))
  expect_identical(dim(estim), c(20L, 7L))
})


test_that("estimate_expectation - glmmTMB", {
  estim <- suppressMessages(estimate_expectation(model2))
  expect_identical(dim(estim), c(644L, 8L))
})


test_that("column name beta regression", {
  data(mtcars)
  mtcars$ord <- datawizard::normalize(mtcars$mpg)
  m <- glmmTMB::glmmTMB(
    ord ~ wt + hp + as.factor(gear) + (1 | cyl),
    data = mtcars,
    family = glmmTMB::ordbeta()
  )
  expect_named(
    estimate_means(m, "gear"),
    c("gear", "Proportion", "SE", "CI_low", "CI_high", "z")
  )
})


test_that("error edge case", {
  set.seed(1)
  d <- data.frame(
    y = rnorm(100),
    x = rnorm(100),
    ID = as.factor(rep(1:50, 2)),
    time = rep(1:2, each = 50),
    int = as.factor(sample(0:1, 100, TRUE)),
    cluster = as.factor(rep(sample.int(8, 50, TRUE), 2))
  )

  m <- glmmTMB::glmmTMB(
    y ~ time * int + x + (1 + time | cluster / ID) + (0 + x | cluster / ID),
    data = d
  )
  expect_error(
    estimate_means(m, "cluster"),
    regex = "You may try using"
  )
})