tests/testthat/test-temporal-frontdoor.R

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

# Temporal RW1 routed through tulpa(): an RW1 random walk is an intrinsic CAR on
# a 1D chain, so it integrates through the nested-Laplace areal (ICAR) path. These
# tests guard the front-door wiring and the fixed-effect recovery.

skip_on_cran()

make_temporal_data <- function(seed, Tt = 24L, m = 25L, b1 = 0.9) {
  set.seed(seed)
  time <- rep(seq_len(Tt), each = m)
  ftrue <- as.numeric(scale(sin(2 * pi * seq_len(Tt) / Tt)))
  x <- rnorm(Tt * m)
  list(df = data.frame(x = x, time = time), time = time, ftrue = ftrue, b1 = b1)
}

# No block in this file reads the outer Pareto-k diagnostic. It runs after the
# grid integration with the RNG state restored, so turning it off leaves the
# fit's coefficients, intervals and evidence unchanged and drops its
# importance-sample inner solves.
.tfd_ctl <- list(diagnose_k = FALSE)

# One binomial RW1 fit through the auto door. Routing, fixed-slope recovery, the
# marginalized SE / interval / prediction and the integrated logLik are separate
# properties of the same fit, so the blocks below read them off this one.
.tfd_rw1_binomial <- fixture_once(function() {
  d <- make_temporal_data(11)
  d$df$y <- rbinom(nrow(d$df), 1, plogis(-0.2 + d$b1 * d$df$x + d$ftrue[d$time]))
  list(d = d, fit = tulpa(y ~ x, data = d$df, family = "binomial",
                          temporal = temporal_rw1("time"), mode = "auto",
                          control = .tfd_ctl))
})

test_that("temporal_rw1 routes through tulpa() and recovers the fixed slope", {
  skip_on_cran()
  fx <- .tfd_rw1_binomial()
  d <- fx$d; fit <- fx$fit

  expect_s3_class(fit, "tulpa_fit")
  expect_equal(fit$backend, "nested_laplace")
  expect_match(fit$selection_reason, "temporal")

  cf <- coef(fit)
  expect_named(cf, c("(Intercept)", "x"))
  expect_lt(abs(cf[["x"]] - d$b1), 0.15)
})

test_that("temporal_rw1 works for a gaussian response", {
  skip_on_cran()
  d <- make_temporal_data(12)
  d$df$y <- -0.2 + d$b1 * d$df$x + d$ftrue[d$time] + rnorm(nrow(d$df), 0, 0.7)
  fit <- tulpa(y ~ x, data = d$df, family = "gaussian",
               temporal = temporal_rw1("time"), mode = "auto", phi = 0.7,
               control = .tfd_ctl)
  expect_lt(abs(coef(fit)[["x"]] - d$b1), 0.15)
})

test_that("nested temporal fit reports a marginalized fixed-effect SE", {
  skip_on_cran()
  fx <- .tfd_rw1_binomial()
  d <- fx$d; fit <- fx$fit

  sm <- coef(summary(fit))
  expect_true(is.finite(sm["x", "std.error"]))   # slope SE from grid marginalization
  ci <- confint(fit)["x", ]
  expect_lt(ci[1], d$b1)
  expect_gt(ci[2], d$b1)

  # se.fit prediction now works (uses the marginalized vcov)
  pr <- predict(fit, newdata = data.frame(x = c(-1, 1)), se.fit = TRUE)
  expect_true(all(is.finite(pr$se.fit)))
})

test_that("nested-fit slope CI covers the truth at roughly nominal rate", {
  skip_if_not_slow()
  hits <- 0L
  for (i in seq_len(15)) {
    d <- make_temporal_data(200 + i)
    d$df$y <- rbinom(nrow(d$df), 1, plogis(-0.2 + d$b1 * d$df$x + d$ftrue[d$time]))
    fit <- tulpa(y ~ x, data = d$df, family = "binomial",
                 temporal = temporal_rw1("time"), mode = "auto",
                 control = .tfd_ctl)
    ci <- confint(fit)["x", ]
    if (ci[1] <= d$b1 && d$b1 <= ci[2]) hits <- hits + 1L
  }
  expect_gte(hits, 12L)
})

test_that("logLik and compare_models return one scalar row per nested fit", {
  skip_on_cran()
  f_t <- .tfd_rw1_binomial()$fit

  expect_length(as.numeric(logLik(f_t)), 1L)   # integrated, not per-grid
  expect_true(is.finite(as.numeric(logLik(f_t))))
  # Two fits reporting the same quantity. A Laplace fit reports a log marginal
  # likelihood and a nested fit a log evidence, and compare_models() refuses to
  # rank those together (test-diagnostics-report-honesty.R).
  cmp <- compare_models(temporal = f_t, again = f_t, criterion = "loglik")
  expect_equal(nrow(cmp), 2L)                   # one row per model, not per grid
})

# RW2 / AR1 route through the native nested-Laplace temporal kernel (issue #73).
test_that("temporal_rw2 routes through tulpa() and recovers the fixed slope", {
  skip_on_cran()
  d <- make_temporal_data(21)
  d$df$y <- rbinom(nrow(d$df), 1, plogis(-0.2 + d$b1 * d$df$x + d$ftrue[d$time]))
  fit <- tulpa(y ~ x, data = d$df, family = "binomial",
               temporal = temporal_rw2("time"), mode = "auto",
               control = .tfd_ctl)
  expect_equal(fit$backend, "nested_laplace")
  expect_match(fit$selection_reason, "rw2")
  expect_lt(abs(coef(fit)[["x"]] - d$b1), 0.15)
})

test_that("temporal_ar1 routes through tulpa() and recovers the fixed slope", {
  skip_on_cran()
  d <- make_temporal_data(31)
  d$df$y <- rbinom(nrow(d$df), 1, plogis(-0.2 + d$b1 * d$df$x + d$ftrue[d$time]))
  fit <- tulpa(y ~ x, data = d$df, family = "binomial",
               temporal = temporal_ar1("time"), mode = "auto",
               control = .tfd_ctl)
  expect_equal(fit$backend, "nested_laplace")
  expect_match(fit$selection_reason, "ar1")
  expect_lt(abs(coef(fit)[["x"]] - d$b1), 0.15)
})

# Additive areal spatial + temporal (space-time) routes through the joint
# multi-block nested-Laplace path (issue #73).
test_that("areal spatial + temporal rw1 routes through the joint nested path", {
  skip_on_cran()
  set.seed(7)
  K <- 12L; Tt <- 10L; m <- 6L
  W <- matrix(0, K, K)
  for (i in seq_len(K - 1L)) { W[i, i + 1L] <- 1; W[i + 1L, i] <- 1 }
  Qs <- diag(rowSums(W)) - W
  eig <- eigen(Qs, symmetric = TRUE)
  phi_s <- as.numeric(eig$vectors[, 1:(K - 1)] %*%
                      (rnorm(K - 1) / sqrt(eig$values[1:(K - 1)])))
  phi_s <- phi_s - mean(phi_s)
  ftrue <- as.numeric(scale(sin(2 * pi * seq_len(Tt) / Tt)))
  grid <- expand.grid(region = seq_len(K), time = seq_len(Tt))
  grid <- grid[rep(seq_len(nrow(grid)), each = m), ]
  grid$x <- rnorm(nrow(grid))
  grid$y <- rbinom(nrow(grid), 1,
                   plogis(-0.2 + 0.9 * grid$x + 0.7 * phi_s[grid$region] +
                          ftrue[grid$time]))
  # The 9 x 9 (tau_spatial x tau_temporal) Cartesian grid warns it is > 50 cells
  # (CCD is a follow-up); not what this test checks.
  fit <- suppressWarnings(
    tulpa(y ~ x + spatial(region), data = grid, family = "binomial",
          spatial = list(type = "icar", adjacency = W),
          temporal = temporal_rw1("time"), mode = "auto"))
  expect_equal(fit$backend, "nested_laplace")
  expect_match(fit$selection_reason, "spatial field \\+ temporal")
  expect_lt(abs(coef(fit)[["x"]] - 0.9), 0.2)
})

# Panel (grouped) temporal: a separate walk per group sharing one tau (issue #73).
test_that("grouped (panel) temporal routes through tulpa()", {
  skip_on_cran()
  set.seed(101)
  G <- 4L; Tt <- 16L; m <- 6L
  shapes <- sapply(seq_len(G), function(g)
    as.numeric(scale(sin(2 * pi * seq_len(Tt) / Tt + g))) * 0.9)
  grid <- expand.grid(time = seq_len(Tt), grp = seq_len(G))
  grid <- grid[rep(seq_len(nrow(grid)), each = m), ]
  grid$grp <- factor(grid$grp); grid$x <- rnorm(nrow(grid))
  grid$y <- rbinom(nrow(grid), 1,
                   plogis(-0.1 + 0.8 * grid$x +
                          shapes[cbind(grid$time, as.integer(grid$grp))]))
  fit <- tulpa(y ~ x, data = grid, family = "binomial",
               temporal = temporal_rw1("time", group_var = "grp"), mode = "auto",
               control = .tfd_ctl)
  expect_s3_class(fit, "tulpa_fit")
  expect_equal(fit$backend, "nested_laplace")
  expect_lt(abs(coef(fit)[["x"]] - 0.8), 0.2)
})

test_that("panel-fit slope CI covers the truth at roughly nominal rate", {
  skip_on_cran()
  b1 <- 0.8; hits <- 0L; nseed <- 12L
  for (i in seq_len(nseed)) {
    set.seed(500 + i)
    G <- 4L; Tt <- 16L; m <- 8L
    shapes <- sapply(seq_len(G), function(g)
      as.numeric(scale(sin(2 * pi * seq_len(Tt) / Tt + g))) * 0.9)
    grid <- expand.grid(time = seq_len(Tt), grp = seq_len(G))
    grid <- grid[rep(seq_len(nrow(grid)), each = m), ]
    grid$grp <- factor(grid$grp); grid$x <- rnorm(nrow(grid))
    grid$y <- rbinom(nrow(grid), 1,
                     plogis(-0.1 + b1 * grid$x +
                            shapes[cbind(grid$time, as.integer(grid$grp))]))
    fit <- tulpa(y ~ x, data = grid, family = "binomial",
                 temporal = temporal_rw1("time", group_var = "grp"), mode = "auto",
                 control = .tfd_ctl)
    ci <- confint(fit)["x", ]
    if (ci[1] <= b1 && b1 <= ci[2]) hits <- hits + 1L
  }
  expect_gte(hits, 9L)   # ~nominal 95% coverage over 12 panel datasets
})

test_that("grouped temporal with one group is identical to the ungrouped fit", {
  skip_on_cran()
  # n_groups == 1 is the degenerate panel case: the flattened node index reduces
  # to the time index and the grouped C++ path must reproduce the single walk
  # bit-for-bit. Guards the panel generalisation against drift.
  set.seed(303)
  Tt <- 18L; m <- 8L
  df <- data.frame(time = rep(seq_len(Tt), each = m), grp = factor(1L))
  ftrue <- as.numeric(scale(sin(2 * pi * seq_len(Tt) / Tt)))
  df$x <- rnorm(nrow(df))
  df$y <- rbinom(nrow(df), 1, plogis(-0.2 + 0.9 * df$x + ftrue[df$time]))
  f_plain   <- tulpa(y ~ x, data = df, family = "binomial",
                     temporal = temporal_rw1("time"), mode = "auto",
                     control = .tfd_ctl)
  f_grouped <- tulpa(y ~ x, data = df, family = "binomial",
                     temporal = temporal_rw1("time", group_var = "grp"),
                     mode = "auto", control = .tfd_ctl)
  expect_equal(as.numeric(logLik(f_grouped)), as.numeric(logLik(f_plain)),
               tolerance = 1e-10)
  expect_equal(unname(coef(f_grouped)), unname(coef(f_plain)), tolerance = 1e-8)
})

test_that("unsupported temporal configurations error clearly", {
  skip_on_cran()
  d <- make_temporal_data(13)
  d$df$y <- rbinom(nrow(d$df), 1, 0.5)

  K <- 6L; W <- matrix(0, K, K)
  for (i in seq_len(K - 1L)) { W[i, i + 1L] <- 1; W[i + 1L, i] <- 1 }
  d$df$region <- factor(sample(seq_len(K), nrow(d$df), replace = TRUE))

  # Panel temporal cannot ride alongside a spatial field (the joint path has no
  # per-group temporal layout yet).
  expect_error(
    tulpa(y ~ x + spatial(region), data = d$df, family = "binomial",
          spatial = list(type = "icar", adjacency = W),
          temporal = temporal_rw1("time", group_var = "region")),
    "panel"
  )
  # A multiscale field has no nested-Laplace kernel, so it cannot host a
  # temporal block through the front door.
  d$df$lon <- rnorm(nrow(d$df)); d$df$lat <- rnorm(nrow(d$df))
  expect_error(
    tulpa(y ~ x, data = d$df, family = "binomial",
          spatial = spatial_multiscale(~ lon + lat),
          temporal = temporal_rw1("time")),
    "areal"
  )
})

test_that("a continuous field shares the nested stack with a temporal field, a smoother and a latent block", {
  skip_on_cran()
  set.seed(8)
  n <- 120L
  d <- data.frame(lon = runif(n), lat = runif(n), time = rep(1:12, 10),
                  x = rnorm(n), z = runif(n))
  d$count <- rpois(n, exp(0.2 + 0.4 * d$x + 0.3 * sin(d$time / 2) +
                            0.4 * sin(3 * d$lon)))
  ctl <- list(axis_refine = "none")
  fit_t <- suppressWarnings(tulpa(
    count ~ x, data = d, family = "poisson", mode = "nested_laplace",
    spatial = spatial_gp(~ lon + lat, nn = 6), temporal = temporal_rw1("time"),
    control = ctl))
  # The same two blocks handed to the registry door directly.
  sp <- tulpa:::validate_gp(spatial_gp(~ lon + lat, nn = 6), d)
  tm <- tulpa:::validate_temporal(temporal_rw1("time"), d)
  direct <- suppressWarnings(tulpa_nested_laplace(
    y = d$count, n_trials = rep(1L, n), X = stats::model.matrix(~ x, d),
    family = "poisson",
    prior = list(tulpa:::.spatial_spec_to_nl_prior(sp),
                 tulpa:::.temporal_spec_to_nl_prior(tm)),
    control = ctl))
  expect_equal(fit_t$backend, "nested_laplace")
  expect_equal(fit_t$log_marginal, direct$log_marginal, tolerance = 1e-8)
  expect_equal(unname(coef(fit_t)), unname(coef(direct)), tolerance = 1e-8)

  fit_s <- suppressWarnings(tulpa(
    count ~ x + s(z), data = d, family = "poisson", mode = "nested_laplace",
    spatial = spatial_gp(~ lon + lat, nn = 6), control = ctl))
  expect_equal(fit_s$backend, "nested_laplace")
  expect_true(all(is.finite(coef(fit_s))))

  blk <- tgmrf(
    Q = function(theta) Matrix::Diagonal(n, exp(-2 * theta[1])),
    prior = function(theta) 0, init = c(log_sigma = 0),
    bounds = list(lower = log(0.2), upper = log(2)), name = "iid")
  fit_l <- suppressWarnings(tulpa(
    count ~ x + latent(blk), data = d, family = "poisson",
    mode = "nested_laplace", spatial = spatial_gp(~ lon + lat, nn = 6),
    control = ctl))
  expect_equal(fit_l$backend, "nested_laplace")
  expect_true(all(is.finite(coef(fit_l))))

  # An HSGP field takes the same stack as a DENSE_BASIS block.
  fit_h <- suppressWarnings(tulpa(
    count ~ x, data = d, family = "poisson", mode = "nested_laplace",
    spatial = spatial_gp(~ lon + lat, approx = "hsgp"),
    temporal = temporal_rw1("time"), control = ctl))
  sh <- tulpa:::validate_hsgp(spatial_gp(~ lon + lat, approx = "hsgp"), d)
  direct_h <- suppressWarnings(tulpa_nested_laplace(
    y = d$count, n_trials = rep(1L, n), X = stats::model.matrix(~ x, d),
    family = "poisson",
    prior = list(tulpa:::.spatial_spec_to_nl_prior(sh),
                 tulpa:::.temporal_spec_to_nl_prior(tm)),
    control = ctl))
  expect_equal(fit_h$backend, "nested_laplace")
  expect_equal(fit_h$log_marginal, direct_h$log_marginal, tolerance = 1e-8)
  expect_equal(unname(coef(fit_h)), unname(coef(direct_h)), tolerance = 1e-8)
})