tests/testthat/test-raster-inference.R

In SSIMmap: The Structural Similarity Index Measure for Maps

test_that("all public transform defaults retain the input values", {
  expect_identical(eval(formals(ssim_raster)$transform)[1], "none")
  expect_identical(eval(formals(ssim_polygon)$transform)[1], "none")
  expect_identical(eval(formals(ssim_bandwidth)$transform)[1], "none")
  a <- terra::rast(nrows = 6, ncols = 6, vals = seq_len(36))
  b <- a * 2 + 5
  expect_identical(ssim_raster(a, b), ssim_raster(a, b, transform = "none"))
  expect_false(isTRUE(all.equal(ssim_raster(a, b)$summary$mean,
                              ssim_raster(a, b, transform = "normal_score")$summary$mean)))
})

# Build the null arrays independently of either inference branch, using only
# the non-testing public calculation for each explicitly shuffled input.
raster_null_reference <- function(a, b, R, seed) {
  av <- terra::values(a, mat = FALSE)
  bv <- terra::values(b, mat = FALSE)
  valid <- which(!is.na(av) & !is.na(bv))
  av[-valid] <- NA_real_
  bv[-valid] <- NA_real_
  terra::values(a) <- av
  terra::values(b) <- bv
  observed <- terra::values(ssim_raster(a, b, global = FALSE))
  null <- array(NA_real_, dim = c(nrow(observed), 4L, R))
  set.seed(seed)
  for (i in seq_len(R)) {
    y <- b
    shuffled <- bv
    shuffled[valid] <- sample(bv[valid], length(valid), replace = FALSE)
    terra::values(y) <- shuffled
    null[, , i] <- terra::values(ssim_raster(a, y, global = FALSE))
  }
  rng <- .Random.seed
  center <- apply(null, c(1, 2), mean)
  local_p <- matrix(NA_real_, nrow(observed), 4L)
  for (j in 1:4) for (cell in which(!is.na(observed[, j]))) {
    z <- null[cell, j, ]
    local_p[cell, j] <- (1 + sum(abs(z - center[cell, j]) >=
      abs(observed[cell, j] - center[cell, j]))) / (R + 1)
  }
  global_null <- t(vapply(seq_len(R), function(i)
    colMeans(null[, , i], na.rm = TRUE), numeric(4)))
  global_observed <- colMeans(observed, na.rm = TRUE)
  global_p <- vapply(1:4, function(j) {
    z <- global_null[, j]
    (1 + sum(abs(z - mean(z)) >= abs(global_observed[j] - mean(z)))) / (R + 1)
  }, numeric(1))
  list(observed = observed, local_p = local_p, global_null = global_null,
       global_p = global_p, rng = rng)
}

test_that("global and local tests use the permutation-centered two-sided null", {
  set.seed(709)
  a <- terra::rast(nrows = 8, ncols = 8, vals = runif(64, 1, 20))
  b <- 21 - a
  terra::values(a)[1] <- NA_real_
  terra::values(b)[64] <- NA_real_
  R <- 39L
  ref <- raster_null_reference(a, b, R, 781)
  g <- ssim_raster(a, b, do_test = TRUE, R = R, seed = 781)
  expect_equal(unname(as.matrix(g$null_distribution)), ref$global_null, tolerance = 1e-12)
  expect_equal(g$summary$p_value, ref$global_p, tolerance = 1e-12)
  expect_equal(g$summary$q_value, p.adjust(ref$global_p, "BH"))
  expect_identical(g$settings$alternative, "two.sided")
  expect_identical(.Random.seed, ref$rng)
  upper_p <- (1 + colSums(sweep(ref$global_null, 2, g$summary$mean, `>=`))) / (R + 1)
  expect_gt(upper_p[4], .9)
  expect_lt(g$summary$p_value[4], .1)

  local <- ssim_raster(a, b, global = FALSE, do_test = TRUE,
                       local_test = TRUE, R = R, seed = 781)
  expect_identical(.Random.seed, ref$rng)
  expect_equal(unname(terra::values(local[[1:4]])), unname(ref$observed))
  expect_equal(unname(terra::values(local[[5:8]])), ref$local_p, tolerance = 1e-12)
  expected_q <- apply(ref$local_p, 2, p.adjust, method = "BH")
  expect_equal(unname(terra::values(local[[9:12]])), expected_q)
  expect_equal(unname(terra::values(local[[13:16]])), (expected_q <= .05) * 1)
  expect_true(all(terra::values(local[[5:12]]) >= 1 / (R + 1), na.rm = TRUE))
  expect_true(all(terra::values(local[[5:12]]) <= 1, na.rm = TRUE))
})

test_that("local replay handles ties, seed=NULL and unadjusted inference", {
  a <- terra::rast(nrows = 6, ncols = 6, vals = rep(4, 36))
  b <- a * 2
  set.seed(19)
  g <- ssim_raster(a, b, do_test = TRUE, R = 7, fdr = FALSE)
  end_state <- .Random.seed
  expect_equal(g$summary$p_value, rep(1, 4))
  set.seed(19)
  local <- ssim_raster(a, b, global = FALSE, do_test = TRUE,
                       local_test = TRUE, R = 7, fdr = FALSE)
  expect_identical(.Random.seed, end_state)
  expect_equal(terra::nlyr(local), 12L)
  expect_true(all(terra::values(local[[5:8]]) == 1, na.rm = TRUE))
  expect_true(all(terra::values(local[[9:12]]) == 0, na.rm = TRUE))
})

test_that("raster significance includes equality at the threshold", {
  a <- terra::rast(nrows = 6, ncols = 6, vals = rep(4, 36))
  b <- a * 2
  for (adjust in c(TRUE, FALSE)) {
    local <- ssim_raster(a, b, global = FALSE, do_test = TRUE,
                         local_test = TRUE, R = 7, fdr = adjust, alpha = 1,
                         seed = 19)
    flags <- local[[paste0(c("SSIM", "SIM", "SIV", "SIP"), "_sig")]]
    expect_true(all(terra::values(flags) == 1, na.rm = TRUE))
  }
})