tests/testthat/test-varcorr.R

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

# Random-effect covariance summary.
#
# The property under test is not just "returns numbers" but that it reports
# WHERE the numbers came from. A Laplace fit conditions on sigma_re -- often the
# default of 1, which the user never chose -- and presenting that beside an
# empirical-Bayes estimate as though they were the same kind of quantity is the
# failure mode this accessor exists to avoid.

vc_data <- function(seed = 1L, G = 30L, per = 12L, sigma = 0.8) {
  set.seed(seed)
  n <- G * per
  g <- rep(seq_len(G), each = per)
  x <- rnorm(n)
  b <- rnorm(G, 0, sigma)
  data.frame(y = rpois(n, exp(0.3 + 0.5 * x + b[g])), x = x, g = factor(g))
}


test_that("a fit with no random effects returns an empty table", {
  set.seed(2)
  d <- data.frame(y = rpois(80, 2), x = rnorm(80))
  fit <- tulpa(y ~ x, data = d, family = "poisson", mode = "laplace")
  vc <- VarCorr(fit)

  expect_s3_class(vc, "data.frame")
  expect_equal(nrow(vc), 0L)
  # Empty but well-formed, so downstream code can bind it without special-casing.
  expect_named(vc, c("term", "coef", "sd", "source"))
})

skip_on_cran_budget()


test_that("a conditioning fit is labelled conditioned, not estimated", {
  skip_on_cran_budget()
  d <- vc_data()
  fit <- tulpa(y ~ x + (1 | g), data = d, family = "poisson",
               mode = "laplace", sigma_re = 0.5)
  vc <- VarCorr(fit)

  expect_equal(vc$source, "conditioned")
  # The value is the input echoed back. That is exactly why it carries a label:
  # without one it is indistinguishable from a fitted 0.5.
  expect_equal(vc$sd, 0.5, tolerance = 1e-8)
})


test_that("the conditioned sigma_re is read off the fit, not the call (#868)", {
  skip_on_cran_budget()
  d <- vc_data()
  # Passed through a variable local to a function: re-evaluating the call
  # outside that frame fell back to 1, or read an unrelated global `s`.
  fitit <- function(dd) {
    s <- 0.3
    tulpa(y ~ x + (1 | g), data = dd, family = "poisson",
          mode = "laplace", sigma_re = s)
  }
  fit <- fitit(d)
  expect_equal(fit$sigma_re_conditioned, 0.3)
  expect_equal(VarCorr(fit)$sd, 0.3, tolerance = 1e-8)
  s <- 5
  expect_equal(VarCorr(fit)$sd, 0.3, tolerance = 1e-8)
  expect_true(any(grepl("0.3", capture.output(print(summary(fit))), fixed = TRUE)))

  # The default it conditioned on is labelled as such.
  fit1 <- suppressWarnings(tulpa(y ~ x + (1 | g), data = d, family = "poisson",
                                 mode = "laplace"))
  expect_equal(VarCorr(fit1)$sd, 1)
  expect_equal(VarCorr(fit1)$source, "conditioned")
})

skip_on_cran()


test_that("an empirical-Bayes fit reports its estimate as estimated", {
  skip_on_cran()
  d <- vc_data()
  fit <- tulpa(y ~ x + (1 | g), data = d, family = "poisson", mode = "eb")
  vc <- VarCorr(fit)

  expect_s3_class(vc, "data.frame")
  expect_equal(nrow(vc), 1L)
  expect_equal(vc$term, "g")
  expect_equal(vc$coef, "(Intercept)")
  expect_equal(vc$source, "estimated")
  # The reported SD is the fit's own estimate, not a re-derivation of it.
  expect_equal(vc$sd, unname(fit$map$sigma), tolerance = 1e-10)
  # And it is in the right neighbourhood of truth, so the wiring is not simply
  # self-consistent about a wrong number.
  expect_lt(abs(vc$sd - 0.8), 0.25)
})


test_that("a sampler fit reports a posterior-mean SD", {
  skip_on_cran()
  d <- vc_data()
  fit <- tulpa(y ~ x + (1 | g), data = d, family = "poisson", mode = "hmc",
               control = list(n_iter = 1200L, warmup = 600L, n_chains = 2L,
                              seed = 4L))
  vc <- VarCorr(fit)

  expect_equal(vc$source, "sampled")
  expect_lt(abs(vc$sd - 0.8), 0.3)

  # Averaged on the SD scale, not the log scale: E[exp(z)] != exp(E[z]), and
  # the log-scale answer is the smaller one for any non-degenerate posterior.
  dm <- tulpa:::.re_draws_mat(fit)
  ls <- dm[, grep("^log_sigma_re", colnames(dm)), drop = FALSE]
  expect_equal(vc$sd, unname(mean(exp(ls[, 1]))), tolerance = 1e-10)
  expect_gt(mean(exp(ls[, 1])), exp(mean(ls[, 1])))
})


test_that("several terms each get a row, in layout order", {
  skip_on_cran()
  set.seed(3)
  G <- 25L; H <- 12L; per <- 14L; n <- G * per
  g <- rep(seq_len(G), each = per)
  h <- sample(seq_len(H), n, replace = TRUE)
  x <- rnorm(n)
  y <- rpois(n, exp(0.3 + 0.4 * x + rnorm(G, 0, 0.9)[g] + rnorm(H, 0, 0.4)[h]))
  d <- data.frame(y = y, x = x, g = factor(g), h = factor(h))

  fit <- tulpa(y ~ x + (1 | g) + (1 | h), data = d, family = "poisson",
               mode = "eb")
  vc <- VarCorr(fit)

  expect_equal(nrow(vc), 2L)
  expect_equal(vc$term, c("g", "h"))
  expect_true(all(vc$sd > 0))
  # The larger simulated component must come back larger; a mis-ordered or
  # recycled assignment would not preserve that.
  expect_gt(vc$sd[1], vc$sd[2])
})


test_that("the covariance matrices ride along with a correlation attribute", {
  skip_on_cran()
  d <- vc_data()
  fit <- tulpa(y ~ x + (1 | g), data = d, family = "poisson", mode = "eb")
  covs <- attr(VarCorr(fit), "cov")

  expect_type(covs, "list")
  expect_named(covs, "g")
  expect_equal(dim(covs$g), c(1L, 1L))
  # A single-coefficient block correlates with itself and nothing else.
  expect_equal(as.numeric(attr(covs$g, "correlation")), 1, tolerance = 1e-10)
  expect_equal(sqrt(as.numeric(covs$g)), VarCorr(fit)$sd, tolerance = 1e-10)
})


test_that("an integrated covariance is reported, not the conditioning fallback", {
  skip_on_cran()
  # The RE-covariance integrators report under `Sigma_mean`, and a fit that
  # integrated Sigma has no sigma_re at all -- so reaching the conditioning
  # fallback invents an input the user never supplied and labels the fit's own
  # result as someone else's assumption (#263). Both backends, because the field
  # is read the same way and only one of them carries `Sigma_draws` alongside.
  set.seed(1)
  site <- rep(seq_len(40), each = 8)
  x    <- rnorm(320)
  b0   <- rnorm(40, 0, 0.7)
  b1   <- rnorm(40, 0, 0.4)
  y    <- rpois(320, exp(0.2 + (0.5 + b1[site]) * x + b0[site]))
  d    <- data.frame(y, x, site = factor(site))

  for (rc in c("gibbs", "nested")) {
    fit <- tulpa(y ~ x + (1 + x | site), data = d, family = "poisson",
                 control = list(re_cov = rc))
    vc  <- VarCorr(fit)

    expect_null(fit$sigma_re, info = rc)
    expect_equal(nrow(vc), 2L, info = rc)
    # re_cov_gibbs SAMPLES Sigma, so its posterior mean is "sampled"; the
    # nested integrator's is "estimated" (gcol33/tulpa#881).
    expect_true(all(vc$source == if (fit$backend == "re_cov_gibbs") "sampled"
                                 else "estimated"), info = rc)
    # The reported SDs are the fit's own integrated covariance, read off
    # Sigma_mean rather than re-derived.
    expect_equal(vc$sd, unname(sqrt(diag(fit$Sigma_mean))),
                 tolerance = 1e-10, info = rc)
    # And near the simulated truth, so the wiring is not self-consistent about
    # a wrong matrix.
    expect_lt(abs(vc$sd[1] - 0.7), 0.25)
    expect_lt(abs(vc$sd[2] - 0.4), 0.25)
  }
})