Add homogeneity_test(): chi2, likelihood-ratio, and Kullback tests for transition regimes

NAMESPACE

@@ -9,6 +9,8 @@ S3method(class_intervals,grd_markov)
 S3method(class_intervals,grd_spatial_markov)
 S3method(first_passage,default)
 S3method(first_passage,grd_markov)
+S3method(homogeneity_test,grd_spatial_markov)
+S3method(homogeneity_test,list)
 S3method(lag_intervals,grd_spatial_markov)
 S3method(mobility_index,default)
 S3method(mobility_index,grd_markov)
@@ -25,6 +27,7 @@ S3method(transition_matrix,grd_spatial_markov)
 export(class_intervals)
 export(classify_dynamics)
 export(first_passage)
+export(homogeneity_test)
 export(lag_intervals)
 export(markov_dynamics)
 export(mobility_index)

NEWS.md

@@ -10,6 +10,18 @@
   `giddy` (`ergodic.mfpt`, `markov.sojourn_time`,
   `mobility.markov_mobility`) via static fixtures.
 
+- `homogeneity_test()` tests whether transition dynamics differ across
+  regimes -- for `grd_spatial_markov` objects, whether the matrices
+  conditioned on spatial-lag class differ from the pooled matrix (Rey, Kang,
+  and Wolf 2016). Reports the Pearson chi-squared Q and likelihood-ratio
+  tests with the Bickenbach--Bald (2003) degrees-of-freedom adjustment and
+  the Kullback (1962) information test. The first two are validated against
+  PySAL `giddy` via static fixtures. The Kullback statistic uses the
+  origin-state margins required by the conditional likelihood, deliberately
+  differing from `giddy.markov.kullback()` through giddy 2.3.9 when transition
+  tables are not flow-balanced. A list method accepts arbitrary regime count
+  matrices.
+
 # griddy 0.1.1
 
 ## Performance

R/inference.R

@@ -0,0 +1,173 @@
+#' Test Homogeneity Of Transition Matrices Across Regimes
+#'
+#' Tests whether transition dynamics are homogeneous across regimes -- for a
+#' spatial Markov object, whether the transition matrices conditioned on
+#' spatial-lag class genuinely differ from the pooled matrix (Rey, Kang, and
+#' Wolf 2016). Three statistics are reported: a Pearson chi-squared Q test and
+#' a likelihood-ratio test with the degrees-of-freedom adjustment of
+#' Bickenbach and Bald (2003), and the Kullback information test (Kullback,
+#' Kupperman, and Ku 1962). The Pearson and likelihood-ratio statistics mirror
+#' the PySAL `giddy` `Spatial_Markov` homogeneity attributes. The Kullback
+#' statistic uses origin-state (row) margins; this deliberately differs from
+#' `giddy.markov.kullback()` through giddy 2.3.9, which uses terminal-state
+#' (column) margins for the stratum-specific term.
+#'
+#' @param x A `grd_spatial_markov` object, or a list of regime transition
+#'   count matrices of identical dimensions (counts, not probabilities).
+#' @param ... Reserved for future methods.
+#'
+#' @return A tibble with one row per test (`chi2`, `likelihood_ratio`,
+#'   `kullback`) and columns `test`, `statistic`, `df`, and `p_value`. The
+#'   null hypothesis is that all regimes share the pooled transition matrix;
+#'   small p-values indicate the conditional dynamics differ.
+#'
+#' @examplesIf identical(Sys.getenv("IN_PKGDOWN"), "true")
+#' panel <- data.frame(
+#'   id = rep(1:9, each = 3),
+#'   year = rep(2020:2022, times = 9),
+#'   value = c(
+#'     1, 2, 3, 4, 3, 5, 6, 7, 6, 8, 9, 10, 11, 10,
+#'     12, 13, 14, 13, 15, 16, 17, 18, 17, 19, 20, 21, 20
+#'   )
+#' )
+#' grid <- sf::st_sf(
+#'   id = 1:9,
+#'   geometry = sf::st_make_grid(
+#'     sf::st_bbox(c(xmin = 0, ymin = 0, xmax = 3, ymax = 3)),
+#'     n = c(3, 3)
+#'   )
+#' ) |>
+#'   dplyr::mutate(
+#'     nb = sfdep::st_contiguity(geometry),
+#'     wt = sfdep::st_weights(nb)
+#'   )
+#' spatial <- spatial_markov(panel, id, year, value, geometry = grid, k = 3)
+#' homogeneity_test(spatial)
+#' @export
+homogeneity_test <- function(x, ...) {
+  UseMethod("homogeneity_test")
+}
+
+#' @export
+homogeneity_test.grd_spatial_markov <- function(x, ...) {
+  mats <- lapply(x$lag_states, function(lag_class) {
+    transition_matrix(x, type = "count", lag_class = lag_class)
+  })
+  names(mats) <- x$lag_states
+  homogeneity_test(mats)
+}
+
+#' @export
+homogeneity_test.list <- function(x, ...) {
+  if (length(x) < 2) {
+    stop("Homogeneity tests need at least two regimes.", call. = FALSE)
+  }
+  mats <- lapply(x, function(m) {
+    m <- as.matrix(m)
+    if (!is.numeric(m)) {
+      stop(
+        "Regime matrices must contain numeric transition counts.",
+        call. = FALSE
+      )
+    }
+    storage.mode(m) <- "double"
+    m
+  })
+  dims <- unique(lapply(mats, dim))
+  if (length(dims) != 1 || dims[[1]][1] != dims[[1]][2]) {
+    stop(
+      "All regime matrices must be square and share the same dimensions.",
+      call. = FALSE
+    )
+  }
+  if (any(vapply(mats, function(m) any(!is.finite(m)), logical(1)))) {
+    stop(
+      "Regime matrices must contain finite transition counts.",
+      call. = FALSE
+    )
+  }
+  if (any(vapply(mats, function(m) any(m < 0), logical(1)))) {
+    stop(
+      "Regime matrices must contain non-negative transition counts.",
+      call. = FALSE
+    )
+  }
+  if (any(vapply(mats, function(m) any(m != floor(m)), logical(1)))) {
+    stop(
+      "Regime matrices must contain whole-number transition counts.",
+      call. = FALSE
+    )
+  }
+  if (any(vapply(mats, function(m) sum(m) == 0, logical(1)))) {
+    stop(
+      "Each regime matrix must contain at least one transition.",
+      call. = FALSE
+    )
+  }
+  bb <- .grd_homogeneity(mats)
+  kb <- .grd_kullback(mats)
+  tibble::tibble(
+    test = c("chi2", "likelihood_ratio", "kullback"),
+    statistic = c(bb$Q, bb$LR, kb$statistic),
+    df = c(bb$dof, bb$dof, kb$dof),
+    p_value = c(bb$Q_p_value, bb$LR_p_value, kb$p_value)
+  )
+}
+
+# Pearson Q and likelihood-ratio homogeneity statistics with the
+# Bickenbach & Bald (2003) degrees-of-freedom adjustment, mirroring
+# giddy.markov.Homogeneity_Results.
+.grd_homogeneity <- function(mats) {
+  r <- nrow(mats[[1]])
+  total <- Reduce(`+`, mats)
+  n_i <- rowSums(total)
+  A_i <- rowSums(total > 0)
+  p_ij <- total / ifelse(n_i == 0, 1, n_i)
+  den <- p_ij + (p_ij == 0)
+
+  b_i <- numeric(r)
+  Q <- 0
+  LR <- 0
+  for (nijm in mats) {
+    nim <- rowSums(nijm)
+    b_i <- b_i + (nim > 0)
+    p_ijm <- nijm / ifelse(nim == 0, 1, nim)
+    Q <- Q + sum(nim * (p_ijm - p_ij)^2 / den)
+    mask <- (nijm > 0) & (p_ij > 0)
+    ratio <- ifelse(mask, p_ijm / p_ij, 1)
+    LR <- LR + sum(nijm * log(ratio))
+  }
+  LR <- 2 * LR
+  dof <- as.integer(sum((b_i - 1) * (A_i - 1)))
+  list(
+    Q = Q,
+    LR = LR,
+    dof = dof,
+    Q_p_value = stats::pchisq(Q, dof, lower.tail = FALSE),
+    LR_p_value = stats::pchisq(LR, dof, lower.tail = FALSE)
+  )
+}
+
+# Kullback information test of conditional homogeneity. The stratum-specific
+# nuisance margins are origin-state (row) totals. giddy <= 2.3.9 instead sums
+# terminal-state columns here; the distinction vanishes for flow-balanced
+# transition tables but not for aggregated panels with endpoint imbalance.
+.grd_kullback <- function(mats) {
+  s <- length(mats)
+  r <- nrow(mats[[1]])
+  xlogx <- function(v) sum(v * log(ifelse(v == 0, 1, v)))
+
+  t1 <- 2 * sum(vapply(mats, xlogx, numeric(1)))
+  pooled <- Reduce(`+`, mats)
+  t2 <- 2 * xlogx(pooled)
+  t3 <- 2 * xlogx(rowSums(pooled))
+  t4 <- 2 * sum(vapply(mats, function(m) xlogx(rowSums(m)), numeric(1)))
+
+  statistic <- max(0, t1 - t4 - t2 + t3)
+  dof <- as.integer(r * (s - 1) * (r - 1))
+  list(
+    statistic = statistic,
+    dof = dof,
+    p_value = stats::pchisq(statistic, dof, lower.tail = FALSE)
+  )
+}

_pkgdown.yml

@@ -54,6 +54,9 @@ reference:
       - sojourn_time
       - first_passage
       - mobility_index
+  - title: Inference
+    contents:
+      - homogeneity_test
   - title: Plotting
     contents:
       - plot_transition_matrix

tests/testthat/fixtures/pysal/pysal_spatial_homogeneity.csv

@@ -0,0 +1,2 @@
+Q,Q_p_value,LR,LR_p_value,dof,kullback,kullback_dof,kullback_p_value
+200.62372763076806,0.0,170.6588284308426,1.492250767398673e-12,60,170.65882843085274,80,1.5963658564297747e-08

tests/testthat/test-inference.R

@@ -0,0 +1,111 @@
+nine_cell_spatial <- function() {
+  panel <- data.frame(
+    id = rep(1:9, each = 3),
+    year = rep(2020:2022, times = 9),
+    value = c(
+      1, 2, 3, 4, 3, 5, 6, 7, 6, 8, 9, 10, 11, 10,
+      12, 13, 14, 13, 15, 16, 17, 18, 17, 19, 20, 21, 20
+    )
+  )
+  grid <- sf::st_sf(
+    id = 1:9,
+    geometry = sf::st_make_grid(
+      sf::st_bbox(c(xmin = 0, ymin = 0, xmax = 3, ymax = 3)),
+      n = c(3, 3)
+    )
+  ) |>
+    dplyr::mutate(
+      nb = sfdep::st_contiguity(geometry),
+      wt = sfdep::st_weights(nb)
+    )
+  spatial_markov(panel, id, year, value, geometry = grid, k = 3)
+}
+
+test_that("kullback statistic reproduces the published Kullback (1962) example", {
+  s1 <- matrix(
+    c(
+      22, 11, 24, 2, 2, 7,
+      5, 23, 15, 3, 42, 6,
+      4, 21, 190, 25, 20, 34,
+      0, 2, 14, 56, 14, 28,
+      32, 15, 20, 10, 56, 14,
+      5, 22, 31, 18, 13, 134
+    ),
+    nrow = 6, byrow = TRUE
+  )
+  s2 <- matrix(
+    c(
+      3, 6, 9, 3, 0, 8,
+      1, 9, 3, 12, 27, 5,
+      2, 9, 208, 32, 5, 18,
+      0, 14, 32, 108, 40, 40,
+      22, 14, 9, 26, 224, 14,
+      1, 5, 13, 53, 13, 116
+    ),
+    nrow = 6, byrow = TRUE
+  )
+
+  out <- homogeneity_test(list(s1, s2))
+  kb <- out[out$test == "kullback", ]
+
+  expect_equal(kb$statistic, 160.961, tolerance = 1e-5)
+  expect_equal(kb$df, 30L)
+  expect_lt(kb$p_value, 1e-6)
+})
+
+test_that("identical regimes are homogeneous under every test", {
+  m <- matrix(c(10, 5, 3, 12), nrow = 2)
+
+  for (scale in c(1, 10, 100, 1000)) {
+    out <- homogeneity_test(list(m * scale, m * scale))
+
+    expect_equal(out$statistic, rep(0, 3), tolerance = 1e-9)
+    expect_equal(out$p_value, rep(1, 3))
+  }
+})
+
+test_that("homogeneity_test returns a tidy three-row tibble", {
+  spatial <- nine_cell_spatial()
+  out <- homogeneity_test(spatial)
+
+  expect_s3_class(out, "tbl_df")
+  expect_equal(out$test, c("chi2", "likelihood_ratio", "kullback"))
+  expect_equal(names(out), c("test", "statistic", "df", "p_value"))
+  expect_true(all(out$statistic >= 0))
+  expect_true(all(out$df >= 0))
+  expect_true(all(out$p_value >= 0 & out$p_value <= 1))
+})
+
+test_that("invalid inputs are rejected", {
+  m <- matrix(c(10, 5, 3, 12), nrow = 2)
+
+  expect_error(homogeneity_test(list(m)), "at least two regimes")
+  expect_error(
+    homogeneity_test(list(m, matrix(0, 3, 3))),
+    "same dimensions"
+  )
+  expect_error(
+    homogeneity_test(list(m, matrix(c(-1, 2, 3, 4), 2))),
+    "non-negative"
+  )
+  expect_error(
+    homogeneity_test(list(m, matrix(c(1, 2, 3, Inf), 2))),
+    "finite"
+  )
+  expect_error(
+    homogeneity_test(list(m, matrix(c(1, 2, 3, NA), 2))),
+    "finite"
+  )
+  expect_error(
+    homogeneity_test(list(m, matrix(c(1, 2, 3, 4.5), 2))),
+    "whole-number"
+  )
+  expect_error(
+    homogeneity_test(list(m, matrix(0, 2, 2))),
+    "at least one transition"
+  )
+  expect_error(
+    homogeneity_test(list(m, matrix(letters[1:4], 2))),
+    "numeric"
+  )
+})