tests/testthat/test-frontdoor-validation.R

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

test_that("tulpa_em_laplace hard-errors on the ignored spatial/re_list (#156g)", {
  e_step <- function(fits, ...) list(weights = 1)
  m_step <- function(weights, ...) list()
  expect_error(
    tulpa_em_laplace(e_step, m_step, spatial = list(type = "icar")),
    "not consumed by tulpa_em_laplace"
  )
  expect_error(
    tulpa_em_laplace(e_step, m_step, re_list = list(list(idx = 1))),
    "not consumed by tulpa_em_laplace"
  )
})

test_that("plot.tulpa_st_summary hard-errors on an unknown type (#156g)", {
  st <- structure(
    data.frame(s = rep(1:3, each = 2), t = rep(1:2, 3), mean = rnorm(6)),
    n_spatial = 3L, n_times = 2L,
    class = c("tulpa_st_summary", "data.frame"))
  expect_error(plot(st, type = "not_a_type"), "Unknown plot")
})

# ============================================================================
# One spatial gate, whichever way the spec arrives (gcol33/tulpa#670, #671)
# ============================================================================

.sp_frontdoor_data <- function(n_units = 10L, n_obs = 100L) {
  set.seed(4L)
  data.frame(s = factor(rep(seq_len(n_units), length.out = n_obs)),
             x = rnorm(n_obs), yb = rbinom(n_obs, 1L, 0.5))
}

.sp_chain_W <- function(n) {
  W <- matrix(0, n, n)
  for (i in seq_len(n - 1L)) { W[i, i + 1L] <- 1; W[i + 1L, i] <- 1 }
  W
}

skip_on_cran_budget()

test_that("a no-intercept single slope (0 + x | g) is rejected on intercept-only paths (#195)", {
  skip_on_cran_budget()
  set.seed(1)
  d <- data.frame(y = rbinom(120, 1, 0.5), x = rnorm(120),
                  g = factor(rep(1:12, 10)))
  # (0 + x | g) has n_coefs == 1 but has_intercept == FALSE: it rides its slope
  # column as a design. The group-index-only paths (AGQ, Gibbs, SPDE, the
  # nested-Laplace native RE) cannot represent that and must reject it clearly
  # rather than silently fit it as a random intercept.
  expect_error(
    tulpa(y ~ x + (0 + x | g), data = d, family = "binomial", mode = "agq"),
    "random-intercept"
  )
  expect_error(
    tulpa(y ~ x + (0 + x | g), data = d, family = "binomial", mode = "gibbs"),
    "random-intercept"
  )
})

test_that("spatiotemporal summary aligns s/t labels with the array when S != T (#197)", {
  skip_on_cran_budget()
  S <- 3L; T <- 2L; nd <- 4L
  # Draws are stored s-fastest: the column for (s, t) is (t - 1) * S + s. Give
  # each column a distinct constant so a mispaired label is visible.
  draws <- matrix(rep(seq_len(S * T), each = nd), nrow = nd)
  obj <- structure(
    list(spatiotemporal = list(n_spatial = S, n_times = T),
         .internal = list(spatiotemporal_draws = draws)),
    class = "tulpa_fit")
  arr  <- spatiotemporal_effects(obj, format = "array")
  summ <- spatiotemporal_effects(obj, format = "summary")
  for (i in seq_len(S)) for (j in seq_len(T)) {
    row <- summ[summ$s == i & summ$t == j, ]
    expect_equal(row$mean, mean(arr[i, j, ]))
  }
})

test_that("ST and latent-factor accessors read a model package's own fit class", {
  skip_on_cran_budget()
  S <- 2L; n_t <- 3L; nd <- 5L
  st_draws <- matrix(rep(seq_len(S * n_t), each = nd), nrow = nd)
  cols <- sprintf("factor[%d,%d]", rep(1:2, 2), rep(1:2, each = 2))
  fit <- structure(
    list(draws = matrix(rnorm(nd * 4), nd, 4, dimnames = list(NULL, cols)),
         spatiotemporal = list(n_spatial = S, n_times = n_t),
         .internal = list(spatiotemporal_draws = st_draws,
                          latent_info = list(type = "latent", n_factors = 2L,
                                             n_obs = 2L))),
    class = "other_pkg_fit")
  summ <- spatiotemporal_effects(fit, format = "summary")
  expect_s3_class(summ, "tulpa_st_summary")
  expect_equal(summ$mean[summ$s == 2 & summ$t == 3], 6)
  lf <- latent_factors(fit)
  expect_equal(nrow(lf), 4L)
  expect_identical(dim(latent_factors(fit, summary = FALSE)), c(nd, 4L))
  expect_error(latent_factors(1:3), "latent factors")
  expect_error(spatiotemporal_effects(1:3), "spatiotemporal")
})

# gcol33/tulpa#613. `.assert_finite_model_inputs()` was wired into three doors by
# hand and not into `tulpa_nested_laplace()`, so an NA response was absorbed
# rather than refused: every grid cell came back `converged = FALSE` with
# `log_marginal = 0`, the softmax over those weights was uniform, a fit was
# returned, and the grid-Hessian retention then reported the user's own missing
# data to them as "This is a tulpa bug; please report." The guard now sits in
# `.validate_glm_design()`, which is what every `(y, X, n_trials)` door goes
# through, so these are one fix rather than four.
test_that("tulpa_nested_laplace() refuses a non-finite response (#613)", {
  skip_on_cran_budget()
  set.seed(1)
  N <- 24L
  region <- rep(seq_len(6L), each = 4L)
  X <- cbind(1, rnorm(N))
  prior <- list(list(type = "iid", obs_idx = region, n_units = 6L,
                     sigma_grid = exp(seq(log(0.2), log(1.5), length.out = 7))))
  call_nl <- function(y, X) {
    tulpa_nested_laplace(
      y = y, n_trials = rep(1L, N), X = X, prior = prior,
      family = "gaussian", phi = 0.49,
      control = list(keep_grid_hessians = TRUE, progress = FALSE))
  }

  # The issue's own repro: an all-NA response.
  expect_error(call_nl(rep(NA_real_, N), X),
               "Non-finite value\\(s\\) in the response")
  # It errors rather than warning about a missing Q, which is what sent the
  # reporter looking for a package defect.
  expect_error(call_nl(rep(NA_real_, N), X), "tulpa_nested_laplace")

  # One bad row, and the message names it.
  y1 <- rnorm(N); y1[7L] <- NA_real_
  expect_error(call_nl(y1, X), "first at row 7")
  y2 <- rnorm(N); y2[3L] <- Inf
  expect_error(call_nl(y2, X), "first at row 3")

  # The model-matrix arm too.
  Xb <- X; Xb[11L, 2L] <- NaN
  expect_error(call_nl(rnorm(N), Xb),
               "Non-finite value\\(s\\) in the model matrix")
  expect_error(call_nl(rnorm(N), Xb), "first at row 11")

  # A clean design still passes the guard (fit cost is tier 2; the guard is
  # what this asserts, so stop at the validator).
  expect_silent(.validate_glm_design(rnorm(N), X, rep(1L, N), "where"))
})

test_that("the guard reaches every door sharing .validate_glm_design (#613)", {
  skip_on_cran_budget()
  set.seed(2)
  N <- 20L
  X <- cbind(1, rnorm(N))
  y <- rbinom(N, 1L, 0.5); y[5L] <- NA_integer_
  expect_error(
    tulpa_gibbs(y = y, n_trials = rep(1L, N), X = X,
                group = rep(1L, N), n_groups = 1L, family = "binomial"),
    "Non-finite value\\(s\\) in the response")
  # The message is prefixed with the door that raised it, not with tulpa().
  expect_error(
    tulpa_gibbs(y = y, n_trials = rep(1L, N), X = X,
                group = rep(1L, N), n_groups = 1L, family = "binomial"),
    "tulpa_gibbs")
  # And the shared validator is where it lives, so a door added later inherits
  # it by construction.
  expect_error(.validate_glm_design(y, X, rep(1L, N), "a_new_door"),
               "a_new_door: Non-finite")
})

test_that("a bare-list adjacency passes the same graph check as a constructor", {
  skip_on_cran_budget()
  d <- .sp_frontdoor_data()
  W <- .sp_chain_W(10L)
  asym <- W
  asym[lower.tri(asym)] <- 0
  expect_error(
    tulpa(yb ~ x + spatial(s), data = d, family = "binomial",
          spatial = list(type = "icar", adjacency = asym), mode = "laplace"),
    "symmetric")
  expect_no_error(
    tulpa(yb ~ x + spatial(s), data = d, family = "binomial",
          spatial = list(type = "icar", adjacency = W), mode = "laplace"))
})

test_that("spatial$type is normalised once, before the selector reads it", {
  skip_on_cran_budget()
  d <- .sp_frontdoor_data()
  W <- .sp_chain_W(10L)
  up <- tulpa(yb ~ x + spatial(s), data = d, family = "binomial",
              spatial = list(type = "ICAR", adjacency = W), mode = "laplace")
  lo <- tulpa(yb ~ x + spatial(s), data = d, family = "binomial",
              spatial = list(type = "icar", adjacency = W), mode = "laplace")
  expect_equal(coef(up), coef(lo))
})

skip_on_cran()
# Front-door input validation (#156f) and hard-errors on silently-ignored
# inputs (#156g). These are cheap structural checks (no fitting), so they run
# on every profile.

test_that("flagship drivers reject a mismatched design (#156f)", {
  skip_on_cran()
  y <- rbinom(20, 1, 0.5)
  region <- rep(1:5, length.out = 20)
  prior <- list(list(type = "iid", obs_idx = region, n_units = 5L,
                     sigma_grid = c(0.5, 1, 2)))
  X <- cbind(1, rnorm(19))                       # nrow(X) != length(y)
  expect_error(
    tulpa_nested_laplace(y = y, n_trials = rep(1L, 20), X = X, prior = prior),
    "nrow\\(X\\)"
  )
  # n_trials length mismatch.
  X2 <- cbind(1, rnorm(20))
  expect_error(
    tulpa_gibbs(y = y, n_trials = rep(1L, 5), X = X2,
                group = rep(1L, 20), n_groups = 1L, family = "binomial"),
    "length\\(n_trials\\)"
  )
  # group out of range.
  expect_error(
    tulpa_gibbs(y = y, n_trials = rep(1L, 20), X = X2,
                group = c(3L, rep(1L, 19)), n_groups = 1L, family = "binomial"),
    "group"
  )
})

test_that("statistical hyperpriors ride re_prior, not control (#156b)", {
  skip_on_cran()
  set.seed(1)
  d <- data.frame(y = rbinom(120, 1, 0.5), x = rnorm(120),
                  g = factor(rep(1:12, 10)))
  # A statistical hyperprior in control is now rejected...
  expect_error(
    tulpa(y ~ x + (1 + x | g), data = d, family = "binomial",
          mode = "laplace", control = list(prior_sigma = c(2, 0.05))),
    "Unknown control knob"
  )
  # ...and a misspelled re_prior key is caught rather than silently ignored.
  expect_error(
    tulpa(y ~ x + (1 | g), data = d, family = "binomial",
          re_prior = list(pri_sigma = 1)),
    "Unknown control knob"
  )
  # A valid re_prior is accepted on the re-covariance path.
  fit <- tulpa(y ~ x + (1 + x | g), data = d, family = "binomial",
               mode = "laplace", re_prior = list(prior_sigma = c(2, 0.05),
                                                  eta = 3))
  expect_identical(fit$backend, "re_cov_nested")
})

test_that("plain Laplace still fits a single random slope (0 + x | g) (#195)", {
  skip_on_cran()
  set.seed(1)
  d <- data.frame(y = rbinom(120, 1, 0.5), x = rnorm(120),
                  g = factor(rep(1:12, 10)))
  # The scalar-sigma_re Laplace path carries the slope column as Z, so it is a
  # valid 1-variance-component fit and must NOT be rejected or redirected to the
  # RE-covariance integrator.
  fit <- tulpa(y ~ x + (0 + x | g), data = d, family = "binomial",
               mode = "laplace")
  expect_identical(fit$backend, "laplace")
})

test_that("the graph is reported once, not once per gate", {
  skip_on_cran()
  d <- .sp_frontdoor_data()
  seen <- character(0)
  withCallingHandlers(
    tulpa(yb ~ x + spatial(s), data = d, family = "binomial",
          spatial = spatial_car(adjacency = .sp_chain_W(10L) * 2, level = "obs"),
          mode = "laplace"),
    warning = function(w) {
      seen <<- c(seen, conditionMessage(w))
      invokeRestart("muffleWarning")
    })
  expect_equal(sum(grepl("weighted graph", seen)), 1L)
})