tests/testthat/test-tulpa-entry.R

In tulpa: Templated Unified Library for Posterior Approximation in Bayesian Hierarchical Models

# test-tulpa-entry.R
# Integration contract for the unified tulpa() entry: formula -> bundle ->
# backend selection -> per-backend arg assembly -> dispatch. Verifies routing,
# the fail-loud reachability contract, and the documented coverage boundaries.

make_pois_re <- function(seed = 1, ng = 8L, per = 10L) {
  set.seed(seed)
  g <- rep(seq_len(ng), each = per)
  x <- rnorm(ng * per)
  eta <- 0.3 + 0.5 * x + rnorm(ng, 0, 0.6)[g]
  data.frame(y = rpois(ng * per, exp(eta)), x = x, g = factor(g))
}

make_binom_re <- function(seed = 3, ng = 8L, per = 12L, nt = 10L) {
  set.seed(seed)
  g <- rep(seq_len(ng), each = per)
  x <- rnorm(ng * per)
  eta <- -0.2 + 0.6 * x + rnorm(ng, 0, 0.5)[g]
  list(
    d = data.frame(y = rbinom(ng * per, nt, plogis(eta)), x = x, g = factor(g)),
    nt = rep(nt, ng * per)
  )
}

skip_on_cran_budget()

test_that("Gibbs rejects unsupported families and non-intercept RE structure", {
  skip_on_cran_budget()
  s <- make_binom_re()
  # gaussian not supported by the PG sampler
  err1 <- expect_error(
    tulpa(y ~ x + (1 | g), s$d, family = "gaussian", mode = "gibbs")
  )
  expect_match(conditionMessage(err1), "binomial.*neg_binomial_2")
  # random slope: more than a single intercept term
  err2 <- expect_error(
    tulpa(y ~ x + (1 + x | g), s$d, family = "binomial", mode = "gibbs",
          n_trials = s$nt)
  )
  expect_match(conditionMessage(err2), "one random-intercept")
})

skip_on_cran()

test_that("tulpa() routes a random-intercept GLMM through the Laplace path", {
  skip_on_cran()
  d <- make_pois_re()
  fit <- tulpa(y ~ x + (1 | g), d, family = "poisson",
               mode = "laplace", sigma_re = 0.6)
  expect_s3_class(fit, "tulpa_fit")
  expect_equal(fit$backend, "laplace")
  expect_equal(fit$inference_tier, 2L)
  expect_identical(fit$family, "poisson")
  expect_true(inherits(fit$formula, "formula"))
})

test_that("tulpa() routes through a sampler (logpost) backend", {
  skip_if_not_slow()
  d <- make_pois_re()
  fit <- tulpa(y ~ x + (1 | g), d, family = "poisson", mode = "mala",
               sigma_re = 0.6, control = list(n_iter = 400L, warmup = 200L))
  expect_s3_class(fit, "tulpa_fit")
  expect_equal(fit$backend, "mala")
  expect_equal(fit$inference_tier, 1L)
})

test_that("tulpa(mode = 'ess') samples a random-effect model jointly (gcol33/tulpa#75)", {
  skip_if_not_slow()
  # The ModelData sampler kernels (ess/hmc/sghmc/...) pack random effects into
  # the latent layout and sample the variance component jointly with the latent
  # and fixed effects -- they are not fixed-effect-only GLMs. ESS only declines
  # the structured spatial / temporal precision (its isotropic proposal cannot
  # carry it; see test-sample-glmm-structure.R); a plain (1 | g) term is sampled.
  d <- make_pois_re()
  fit <- suppressWarnings(
    tulpa(y ~ x + (1 | g), d, family = "poisson", mode = "ess",
          control = list(n_iter = 300L, warmup = 150L)))
  expect_s3_class(fit, "tulpa_fit")
  expect_equal(fit$backend, "ess")
  expect_equal(fit$inference_tier, 1L)
  # The RE variance is sampled (full Bayes), not conditioned away.
  expect_true("log_sigma_re" %in% fit$param_names)
  expect_true(all(is.finite(fit$means)))
})

test_that("tulpa(mode = 'ess'/'hmc') fits a fixed-effect GLM end to end", {
  skip_on_cran()
  set.seed(3)
  n <- 400L; x <- rnorm(n)
  d <- data.frame(y = rbinom(n, 1L, plogis(-0.3 + 0.8 * x)), x = x)
  for (m in c("ess", "hmc")) {
    ctl <- list(n_iter = 1000L, warmup = 500L)
    # n_chains is read by the NUTS kernel only; every other backend returns one
    # set of draws and now refuses it rather than ignoring it
    # (gcol33/tulpa#704).
    if (m == "hmc") ctl$n_chains <- 2L
    fit <- tulpa(y ~ x, d, family = "binomial", mode = m, control = ctl)
    expect_s3_class(fit, "tulpa_fit")
    expect_equal(fit$backend, m)
    expect_true(all(is.finite(fit$means)))
  }
  expect_error(
    tulpa(y ~ x, d, family = "binomial", mode = "ess",
          control = list(n_iter = 100L, warmup = 50L, n_chains = 2L)),
    "only read by the NUTS")
})

test_that("correlated random slopes: logpost path works, Laplace path integrates Sigma", {
  skip_if_not_slow()
  d <- make_pois_re(seed = 2)
  # A correlated term has no scalar sigma_re to condition on, so the sampler
  # modes redirect to the Metropolis-within-Gibbs debias exactly as the
  # deterministic mode redirects to the nested integrator; `sigma_re` is then
  # unused and warns.
  fit <- suppressWarnings(
    tulpa(y ~ x + (1 + x | g), d, family = "poisson", mode = "mala",
          sigma_re = 0.5, control = list(n_iter = 300L, warmup = 150L)))
  expect_equal(fit$backend, "re_cov_gibbs")
  # Laplace path now routes a single correlated `(1 + x | g)` term to the
  # nested-Laplace Sigma integrator (no scalar sigma_re to condition on).
  fit2 <- tulpa(y ~ x + (1 + x | g), d, family = "poisson", mode = "laplace",
                control = list(seed = 1L, n_draws = 300L))
  expect_equal(fit2$backend, "re_cov_nested")
  expect_equal(fit2$inference_tier, 2L)
  expect_setequal(fit2$posterior$parameter,
                  c("sigma_1", "sigma_2", "rho_12",
                    "Sigma_11", "Sigma_12", "Sigma_22"))
  # the exact debias is reachable via control$re_cov
  fit3 <- tulpa(y ~ x + (1 + x | g), d, family = "poisson", mode = "laplace",
                control = list(re_cov = "gibbs", n_iter = 400L, warmup = 200L,
                               seed = 1L))
  expect_equal(fit3$backend, "re_cov_gibbs")
})

test_that("tulpa() fits a single random slope (0 + x | g) on the Laplace path", {
  skip_on_cran()
  set.seed(91L)
  G <- 40L; npg <- 12L; N <- G * npg
  g <- rep(seq_len(G), each = npg); x <- rnorm(N)
  u <- rnorm(G, 0, 0.7)
  d <- data.frame(y = rpois(N, exp(0.3 + 0.5 * x + x * u[g])), x = x,
                  g = factor(g))
  # n_coefs == 1 (a single random slope, no intercept) stays on the scalar
  # sigma_re design path -- it does not trigger the RE-covariance redirect.
  fit <- suppressMessages(
    tulpa(y ~ x + (0 + x | g), d, family = "poisson", mode = "laplace",
          sigma_re = 0.7))
  expect_s3_class(fit, "tulpa_fit")
  expect_equal(fit$backend, "laplace")
  expect_equal(fit$inference_tier, 2L)
})


test_that("tulpa() validates family and defaults sigma_re with a message", {
  skip_on_cran()
  d <- make_pois_re()
  expect_error(tulpa(y ~ x + (1 | g), d, family = "weibull"), "Unknown family")
  expect_warning(
    tulpa(y ~ x + (1 | g), d, family = "poisson", mode = "laplace"),
    "sigma_re = 1"
  )
})

test_that("auto mode infers a random-intercept's SD instead of conditioning it at 1 (gcol33/tulpa#267)", {
  # A random slope already routes auto to the RE-covariance integrator (see
  # "correlated random slopes" above); a bare `(1 | g)` term used to fall
  # through to the plain-GLM default (backend "mala") and condition sigma_re
  # at 1 instead. auto now treats the two the same way.
  skip_on_cran()
  set.seed(2)
  G <- 15L; npg <- 30L; n <- G * npg
  g <- factor(rep(seq_len(G), each = npg)); x <- rnorm(n)
  b0 <- rnorm(G, 0, 0.9)
  y <- rpois(n, exp(0.2 + 0.3 * x + b0[as.integer(g)]))
  d <- data.frame(y = y, x = x, g = g)

  fit <- tulpa(y ~ x + (1 | g), data = d, family = "poisson",
               control = list(n_iter = 1500L, warmup = 800L, seed = 5L))
  expect_equal(fit$backend, "re_cov_gibbs")
  expect_equal(fit$inference_tier, 1L)

  vc <- VarCorr(fit)
  # Inferred, not conditioned: re_cov_gibbs samples Sigma (gcol33/tulpa#881).
  expect_equal(vc$source, "sampled")
  # A real recovery bound against the simulated truth (0.9), not a shape check.
  expect_lt(abs(vc$sd - 0.9), 0.35)

  # An explicit conditional mode is unaffected: it still conditions (and warns)
  # rather than being silently redirected out from under the caller.
  expect_warning(
    fit_laplace <- tulpa(y ~ x + (1 | g), data = d, family = "poisson",
                         mode = "laplace"),
    "sigma_re = 1"
  )
  expect_equal(fit_laplace$backend, "laplace")
  expect_equal(VarCorr(fit_laplace)$source, "conditioned")
})

test_that("tulpa() handles a no-RE model on the design path", {
  skip_on_cran()
  set.seed(5)
  d <- data.frame(y = rpois(120, exp(0.4 + 0.3 * rnorm(120))), x = rnorm(120))
  fit <- tulpa(y ~ x, d, family = "poisson", mode = "laplace")
  expect_s3_class(fit, "tulpa_fit")
  expect_equal(fit$backend, "laplace")
})

test_that("tulpa() routes through imh_laplace (Laplace proposal + MH)", {
  skip_if_not_slow()
  d <- make_pois_re()
  fit <- tulpa(y ~ x + (1 | g), d, family = "poisson", mode = "imh_laplace",
               sigma_re = 0.6, control = list(n_iter = 400L, warmup = 200L))
  expect_s3_class(fit, "tulpa_fit")
  expect_equal(fit$backend, "imh_laplace")
  expect_equal(fit$inference_tier, 1L)
  expect_true(is.finite(fit$mean_accept))
})

test_that("tulpa() routes a binomial random-intercept model through Gibbs", {
  skip_if_not_slow()
  s <- make_binom_re()
  fit <- tulpa(y ~ x + (1 | g), s$d, family = "binomial", mode = "gibbs",
               n_trials = s$nt, control = list(n_iter = 400L, warmup = 200L))
  expect_s3_class(fit, "tulpa_fit")
  expect_equal(fit$backend, "gibbs")
  expect_equal(fit$inference_tier, 1L)
})