tests/testthat/test-polygon-inference.R

In SSIMmap: The Structural Similarity Index Measure for Maps

polygon_inference_fixture <- function() {
  xy <- expand.grid(x = 0:3, y = 0:2)
  xy$x <- xy$x + seq_len(nrow(xy))^2 / 1000
  geometry <- sf::st_sfc(lapply(seq_len(nrow(xy)), function(i) {
    x <- xy$x[i]; y <- xy$y[i]
    sf::st_polygon(list(matrix(c(x, y, x + .2, y, x + .2, y + .2,
                                x, y + .2, x, y), ncol = 2, byrow = TRUE)))
  }), crs = 3857)
  a <- c(2, 19, 7, 31, 11, 43, 23, 5, 37, 13, 29, 17)
  sf::st_sf(a = a, b = 50 - a, geometry = geometry)
}

# Construct reference null distributions using only the descriptive public
# calculation on explicitly shuffled maps, independently of inference code.
polygon_null_reference <- function(shape, R, seed) {
  metrics <- c("SSIM", "SIM", "SIV", "SIP")
  observed <- as.matrix(sf::st_drop_geometry(ssim_polygon(
    shape, "a", "b", global = FALSE, bandwidth = 4, transform = "none"))[metrics])
  null <- array(NA_real_, c(nrow(shape), 4L, R))
  set.seed(seed)
  for (i in seq_len(R)) {
    shuffled <- shape
    shuffled$a <- sample(shape$a, replace = FALSE)
    shuffled$b <- sample(shape$b, replace = FALSE)
    null[, , i] <- as.matrix(sf::st_drop_geometry(ssim_polygon(
      shuffled, "a", "b", global = FALSE, bandwidth = 4,
      transform = "none"))[metrics])
  }
  rng <- .Random.seed
  p <- matrix(NA_real_, nrow(shape), 4L, dimnames = list(NULL, metrics))
  for (j in seq_along(metrics)) for (i in seq_len(nrow(shape))) {
    z <- null[i, j, ]; mu <- mean(z)
    p[i, j] <- (1 + sum(abs(z - mu) >= abs(observed[i, j] - mu))) / (R + 1)
  }
  global_null <- t(vapply(seq_len(R), function(i) colMeans(null[, , i]), numeric(4)))
  global_observed <- colMeans(observed)
  global_p <- vapply(seq_along(metrics), function(j) {
    z <- global_null[, j]; mu <- mean(z)
    (1 + sum(abs(z - mu) >= abs(global_observed[j] - mu))) / (R + 1)
  }, numeric(1))
  list(observed = observed, p = p, global_null = global_null,
       global_p = global_p, rng = rng)
}

test_that("polygon local components match independent two-sided permutation tests", {
  shape <- polygon_inference_fixture()
  original <- shape
  metrics <- c("SSIM", "SIM", "SIV", "SIP")
  ref <- polygon_null_reference(shape, R = 39L, seed = 804)
  result <- ssim_polygon(shape, "a", "b", global = FALSE, bandwidth = 4,
                         transform = "none", do_test = TRUE, R = 39, seed = 804)
  expect_identical(.Random.seed, ref$rng)
  expect_identical(shape, original)
  expect_identical(sf::st_geometry(result), sf::st_geometry(shape))
  values <- sf::st_drop_geometry(result)
  expect_equal(as.matrix(values[metrics]), ref$observed, tolerance = 1e-12)
  expected_q <- apply(ref$p, 2, p.adjust, method = "BH")
  for (m in metrics) {
    expect_equal(values[[paste0(m, "_p")]], ref$p[, m], tolerance = 1e-12)
    expect_equal(values[[paste0(m, "_q")]], expected_q[, m], tolerance = 1e-12)
    expect_identical(values[[paste0(m, "_sig")]], expected_q[, m] <= .05)
  }
  expect_identical(values$p_value, values$SSIM_p)
  expect_identical(values$q_value, values$SSIM_q)
  expect_identical(values$sig, values$SSIM_sig)
  expect_true(all(ref$p >= 1 / 40 & ref$p <= 1))
  expect_gt(length(unique(unlist(values[paste0(metrics, "_p")]))), 1L)

  invisible(capture.output(global <- ssim_polygon(shape, "a", "b",
    global = TRUE, bandwidth = 4, transform = "none", do_test = TRUE,
    R = 39, seed = 804)))
  expect_equal(global$p_global$P_value, ref$global_p, tolerance = 1e-12)
  expect_identical(.Random.seed, ref$rng)
  upper_sip <- (1 + sum(ref$global_null[, 4] >= mean(ref$observed[, 4]))) / 40
  expect_gt(upper_sip, .9)
  expect_lt(global$p_global$P_value[4], .1)
})

test_that("polygon thresholds include equality with and without FDR adjustment", {
  shape <- polygon_inference_fixture()
  metrics <- c("SSIM", "SIM", "SIV", "SIP")
  for (adjust in c(TRUE, FALSE)) {
    initial <- ssim_polygon(shape, "a", "b", global = FALSE, bandwidth = 4,
      do_test = TRUE, R = 19, seed = 804, fdr = adjust)
    threshold <- initial$SIP_q[1]
    expect_gt(threshold, 0)
    expect_lt(threshold, 1)
    at_boundary <- ssim_polygon(shape, "a", "b", global = FALSE, bandwidth = 4,
      do_test = TRUE, R = 19, seed = 804, fdr = adjust, alpha = threshold)
    expect_true(at_boundary$SIP_sig[1])
    for (m in metrics) {
      expect_identical(at_boundary[[paste0(m, "_sig")]],
                       at_boundary[[paste0(m, "_q")]] <= threshold)
      if (!adjust) expect_identical(at_boundary[[paste0(m, "_q")]],
                                    at_boundary[[paste0(m, "_p")]])
    }
  }
})

test_that("polygon ties and seed=NULL retain all four local tests", {
  shape <- polygon_inference_fixture()
  shape$a <- 4; shape$b <- 8
  set.seed(92)
  result <- ssim_polygon(shape, "a", "b", global = FALSE, bandwidth = 4,
                         do_test = TRUE, R = 7, fdr = TRUE, alpha = 1)
  for (m in c("SSIM", "SIM", "SIV", "SIP")) {
    expect_equal(result[[paste0(m, "_p")]], rep(1, nrow(shape)))
    expect_equal(result[[paste0(m, "_q")]], rep(1, nrow(shape)))
    expect_true(all(result[[paste0(m, "_sig")]]))
  }
  plain <- ssim_polygon(shape, "a", "b", global = FALSE, bandwidth = 4)
  expect_identical(names(plain), c(names(shape), "SSIM", "SIM", "SIV", "SIP"))
})