tests/testthat/test-compatibility-and-sampling.R

In MRS: Multi-Resolution Scanning for Cross-Sample and Cross-Group Differences

test_that("positive default bounds retain the 1.2.6 definition", {
  X <- cbind(seq(-2, 3, length.out = 40), seq(0.1, 4, length.out = 40))
  G <- rep(1:2, each = 20)
  expected <- t(apply(X, 2L, range))
  expected[, 2L] <- expected[, 2L] * 1.0001

  fit <- mrs(X, G, K = 2)
  expect_equal(fit$Data$Omega, expected, tolerance = 0)
})

test_that("ANDOVA defaults reproduce a deterministic 1.2.6 result", {
  X <- matrix(c(seq(-1.6, 1.3, length.out = 24),
                seq(-1.1, 1.8, length.out = 24)), ncol = 1L)
  G <- rep(1:2, each = 24L)
  H <- rep(rep(1:3, each = 8L), 2L)

  fit <- andova(X, G, H, K = 2, nu_vec = c(0.5, 1))

  expect_equal(fit$PostGlobNull, 0.763795901666116, tolerance = 1e-12)
  expect_equal(fit$PriorGlobNull, 0.694184400112226, tolerance = 1e-12)
  expect_equal(fit$LogLikelihood, 14.4390440353911, tolerance = 1e-12)
  expect_equal(
    fit$RepresentativeTree$AltProbs,
    c(0.0464317065424557, 0.0469209729815043, 0.0694544545024379,
      0.0469210053986832, 0.0464317182535477, 0.0694544838993035,
      0.126264487972054),
    tolerance = 1e-12
  )
  expect_equal(fit$Data$Omega, matrix(c(-1.6, 1.80018), nrow = 1L),
               tolerance = 1e-15)
})

test_that("multidimensional MRS posterior trees are valid and reproducible", {
  set.seed(12345)
  X <- matrix(rnorm(180), ncol = 3L)
  G <- rep(1:3, each = 20L)

  set.seed(12345)
  fit <- mrs(X, G, K = 2, n_post_samples = 25)
  set.seed(12345)
  repeat_fit <- mrs(X, G, K = 2, n_post_samples = 25)

  expect_equal(fit$PostSamples, repeat_fit$PostSamples)
  expect_length(fit$PostSamples, 25L)
  expect_true(all(vapply(fit$PostSamples, function(tree) {
    nrow(tree$Regions) == length(tree$Levels) &&
      ncol(tree$Regions) == 2L * ncol(X) &&
      nrow(tree$EffectSizes) == length(tree$Levels) &&
      ncol(tree$EffectSizes) == length(unique(G)) &&
      all(is.finite(tree$EffectSizes)) &&
      all(tree$Directions %in% seq_len(ncol(X)))
  }, logical(1L))))
  expect_setequal(unique(unlist(lapply(fit$PostSamples, `[[`, "Directions"))),
                  seq_len(ncol(X)))
})

test_that("ANDOVA posterior sampling stays inside the probability domain", {
  set.seed(12345)
  X <- matrix(rnorm(144), ncol = 2L)
  G <- rep(1:3, each = 24L)
  H <- rep(rep(1:3, each = 8L), 3L)

  for (method in c("newton", "riemann")) {
    set.seed(12345)
    fit <- andova(X, G, H, K = 2, nu_vec = c(0.2, 1, 5),
                  n_post_samples = 100, method = method, n_grid_theta = 20)
    effect_sizes <- unlist(lapply(fit$PostSamples, `[[`, "EffectSizes"))
    expect_length(fit$PostSamples, 100L)
    expect_true(all(is.finite(effect_sizes)), info = method)
    expect_true(all(unlist(lapply(fit$PostSamples, `[[`, "Directions")) %in% 1:2),
                info = method)
  }
})

test_that("posterior samples do not require returning a representative tree", {
  set.seed(12345)
  X <- matrix(rnorm(80), ncol = 2L)
  G <- rep(1:2, each = 20L)
  fit <- mrs(X, G, K = 2, return_tree = FALSE, n_post_samples = 3)

  expect_length(fit$RepresentativeTree, 0L)
  expect_length(fit$PostSamples, 3L)
  expect_true(all(vapply(fit$PostSamples, function(tree) {
    all(is.finite(tree$EffectSizes))
  }, logical(1L))))
})

test_that("fit probabilities, regions, and observation indices obey invariants", {
  set.seed(12345)
  X <- rbind(matrix(rnorm(120), ncol = 2L),
             matrix(rnorm(120, 0.4), ncol = 2L))
  G <- rep(1:2, each = 60L)
  fit <- mrs(X, G, K = 4)
  tree <- fit$RepresentativeTree

  expect_true(fit$PostGlobNull >= 0 && fit$PostGlobNull <= 1)
  expect_true(fit$PriorGlobNull >= 0 && fit$PriorGlobNull <= 1)
  expect_true(all(tree$AltProbs >= 0 & tree$AltProbs <= 1))
  expect_true(all(is.finite(tree$EffectSizes)))
  expect_true(all(tree$Levels >= 0 & tree$Levels <= fit$OtherInfo$K))
  expect_true(all(tree$Directions %in% seq_len(ncol(X))))
  expect_true(all(vapply(tree$DataPoints, function(index) {
    all(index >= 1 & index <= nrow(X)) && !anyDuplicated(index)
  }, logical(1L))))

  for (dimension in seq_len(ncol(X))) {
    lower <- tree$Regions[, 2L * dimension - 1L]
    upper <- tree$Regions[, 2L * dimension]
    expect_true(all(lower < upper))
    expect_true(all(lower >= fit$Data$Omega[dimension, 1L]))
    expect_true(all(upper <= fit$Data$Omega[dimension, 2L]))
  }
})