From 0f6bb96868ae2e10a20b028dd5eba2982ab890be Mon Sep 17 00:00:00 2001
From: Joe Thorley <joe@poissonconsulting.ca>
Date: Tue, 25 Aug 2026 10:19:11 +0100
Subject: [PATCH 1/3] Test the multinomial functions against their binomial
 special case

A two-category multinomial is a binomial, which gives an independent
check on all four functions: trial-summed log-likelihoods match
log_lik_binom(), trial-summed deviances (and summed squared deviance
residuals) match dev_binom(), the standardized residual of the first
category matches res_binom() with the second its negative, and a
two-category draw is stream-identical to ran_binom() with binomial
marginal mean and variance.

Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
---
 tests/testthat/test-dev.R     | 26 ++++++++++++++++++++
 tests/testthat/test-log-lik.R | 17 +++++++++++++
 tests/testthat/test-ran.R     | 31 +++++++++++++++++++++++
 tests/testthat/test-res.R     | 46 +++++++++++++++++++++++++++++++++++
 4 files changed, 120 insertions(+)

diff --git a/tests/testthat/test-dev.R b/tests/testthat/test-dev.R
index b4fcb0cc..5b987712 100644
--- a/tests/testthat/test-dev.R
+++ b/tests/testthat/test-dev.R
@@ -569,6 +569,32 @@ test_that("dev_multinom res", {
   expect_equal(sign(res), sign(x - size * prob))
   expect_equal(sum(res^2), sum(dev_multinom(x, size, prob, group)))
 })
+test_that("dev_multinom with two categories matches dev_binom", {
+  x <- c(0, 3, 7, 10)
+  size <- 10
+  prob <- c(0.05, 0.2, 0.5, 0.9)
+  group <- rep(seq_along(x), each = 2)
+  x_long <- as.vector(rbind(x, size - x))
+  prob_long <- as.vector(rbind(prob, 1 - prob))
+
+  # the Poisson-form rows carry offsets that the binomial deviance doesn't,
+  # but they cancel within a trial, so the trial totals agree exactly
+  dev <- dev_multinom(x_long, size = size, prob = prob_long, group = group)
+  expect_equal(colSums(matrix(dev, nrow = 2)), dev_binom(x, size, prob))
+
+  res <- dev_multinom(
+    x_long,
+    size = size,
+    prob = prob_long,
+    group = group,
+    res = TRUE
+  )
+  expect_equal(colSums(matrix(res, nrow = 2)^2), dev_binom(x, size, prob))
+  expect_equal(
+    sign(matrix(res, nrow = 2)[1, ]),
+    sign(res_binom(x, size, prob, type = "dev"))
+  )
+})
 
 test_that("dev_neg_binom", {
   expect_identical(
diff --git a/tests/testthat/test-log-lik.R b/tests/testthat/test-log-lik.R
index 37abfd1d..46518591 100644
--- a/tests/testthat/test-log-lik.R
+++ b/tests/testthat/test-log-lik.R
@@ -427,6 +427,23 @@ test_that("log_lik_multinom", {
     dmultinom(x[group == 2], size = 4, prob = prob[group == 2], log = TRUE)
   )
 })
+test_that("log_lik_multinom with two categories matches log_lik_binom", {
+  # a two-category multinomial is a binomial, so a trial's log-likelihood
+  # (the sum over its rows) must match log_lik_binom() on the first category
+  x <- c(0, 3, 7, 10)
+  size <- 10
+  prob <- c(0.05, 0.2, 0.5, 0.9)
+  log_lik <- log_lik_multinom(
+    as.vector(rbind(x, size - x)),
+    size = size,
+    prob = as.vector(rbind(prob, 1 - prob)),
+    group = rep(seq_along(x), each = 2)
+  )
+  expect_equal(
+    colSums(matrix(log_lik, nrow = 2)),
+    log_lik_binom(x, size, prob)
+  )
+})
 
 test_that("log_lik_neg_binom", {
   expect_identical(
diff --git a/tests/testthat/test-ran.R b/tests/testthat/test-ran.R
index 68b11318..4b9886d6 100644
--- a/tests/testthat/test-ran.R
+++ b/tests/testthat/test-ran.R
@@ -196,6 +196,37 @@ test_that("ran_multinom", {
     expect_identical(x[3] + x[4], 6L)
   })
 })
+test_that("ran_multinom with two categories matches ran_binom", {
+  size <- c(10, 4, 20)
+  prob <- c(0.2, 0.5, 0.75)
+  # rmultinom() draws a two-category trial with a single rbinom() call, so
+  # the first category is stream-identical to ran_binom()
+  withr::with_seed(42, {
+    x <- ran_multinom(
+      size = rep(size, each = 2),
+      prob = as.vector(rbind(prob, 1 - prob)),
+      group = rep(seq_along(size), each = 2)
+    )
+  })
+  withr::with_seed(42, {
+    y <- ran_binom(length(size), size = size, prob = prob)
+  })
+  expect_identical(matrix(x, nrow = 2)[1, ], y)
+
+  # and the marginal of a category is binomial, so its mean and variance
+  # match size * prob and size * prob * (1 - prob)
+  n_group <- 20000
+  withr::with_seed(7, {
+    x <- ran_multinom(
+      size = rep(10, 2 * n_group),
+      prob = rep(c(0.3, 0.7), n_group),
+      group = rep(seq_len(n_group), each = 2)
+    )
+  })
+  first <- matrix(x, nrow = 2)[1, ]
+  expect_equal(mean(first), 10 * 0.3, tolerance = 0.01)
+  expect_equal(var(first), 10 * 0.3 * 0.7, tolerance = 0.01)
+})
 
 test_that("ran_neg_binom", {
   expect_error(ran_neg_binom(NA_integer_))
diff --git a/tests/testthat/test-res.R b/tests/testthat/test-res.R
index 0fbf9c5b..4cc11645 100644
--- a/tests/testthat/test-res.R
+++ b/tests/testthat/test-res.R
@@ -649,6 +649,52 @@ test_that("res_multinom simulate", {
     expect_equal(sd(res), 1.00164264126585)
   })
 })
+test_that("res_multinom with two categories matches res_binom", {
+  x <- c(0, 3, 7, 10)
+  size <- 10
+  prob <- c(0.05, 0.2, 0.5, 0.9)
+  group <- rep(seq_along(x), each = 2)
+  x_long <- as.vector(rbind(x, size - x))
+  prob_long <- as.vector(rbind(prob, 1 - prob))
+
+  res_type <- function(type) {
+    matrix(
+      res_multinom(x_long, size, prob_long, group, type = type),
+      nrow = 2
+    )
+  }
+
+  expect_equal(res_type("raw")[1, ], res_binom(x, size, prob, type = "raw"))
+  # the second category is the first one's complement, so its standardized
+  # residual is the negative of it
+  standardized <- res_binom(x, size, prob, type = "standardized")
+  expect_equal(res_type("standardized")[1, ], standardized)
+  expect_equal(res_type("standardized")[2, ], -standardized)
+  # a per-category deviance residual isn't the binomial one, but squaring
+  # and summing within the trial recovers the binomial deviance
+  expect_equal(colSums(res_type("dev")^2), dev_binom(x, size, prob))
+  expect_equal(
+    sign(res_type("dev")[1, ]),
+    sign(res_binom(x, size, prob, type = "dev"))
+  )
+
+  # a simulated two-category trial is a single rbinom() draw, so it is
+  # stream-identical to res_binom(simulate = TRUE)
+  withr::with_seed(3, {
+    sim <- res_multinom(
+      x_long,
+      size,
+      prob_long,
+      group,
+      type = "data",
+      simulate = TRUE
+    )
+  })
+  withr::with_seed(3, {
+    sim_binom <- res_binom(x, size, prob, type = "data", simulate = TRUE)
+  })
+  expect_equal(matrix(sim, nrow = 2)[1, ], sim_binom)
+})
 
 test_that("res_neg_binom", {
   expect_identical(

From 155ab3143aab0e30e0051972c2b6c87bbb6628b5 Mon Sep 17 00:00:00 2001
From: Joe Thorley <joe@poissonconsulting.ca>
Date: Tue, 25 Aug 2026 10:22:55 +0100
Subject: [PATCH 2/3] Drop the RNG-stream assertions from the binomial checks

Stream-identity to rbinom() relies on rmultinom() internals rather than a
documented guarantee. ran_multinom() now checks the first category's mean,
variance and count distribution against the binomial instead, and the
res_multinom() simulate check is dropped as redundant with it.

Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
---
 tests/testthat/test-dev.R     |  1 +
 tests/testthat/test-log-lik.R |  1 +
 tests/testthat/test-ran.R     | 33 +++++++++++----------------------
 tests/testthat/test-res.R     | 18 +-----------------
 4 files changed, 14 insertions(+), 39 deletions(-)

diff --git a/tests/testthat/test-dev.R b/tests/testthat/test-dev.R
index 5b987712..42e597de 100644
--- a/tests/testthat/test-dev.R
+++ b/tests/testthat/test-dev.R
@@ -569,6 +569,7 @@ test_that("dev_multinom res", {
   expect_equal(sign(res), sign(x - size * prob))
   expect_equal(sum(res^2), sum(dev_multinom(x, size, prob, group)))
 })
+
 test_that("dev_multinom with two categories matches dev_binom", {
   x <- c(0, 3, 7, 10)
   size <- 10
diff --git a/tests/testthat/test-log-lik.R b/tests/testthat/test-log-lik.R
index 46518591..3d808ef0 100644
--- a/tests/testthat/test-log-lik.R
+++ b/tests/testthat/test-log-lik.R
@@ -427,6 +427,7 @@ test_that("log_lik_multinom", {
     dmultinom(x[group == 2], size = 4, prob = prob[group == 2], log = TRUE)
   )
 })
+
 test_that("log_lik_multinom with two categories matches log_lik_binom", {
   # a two-category multinomial is a binomial, so a trial's log-likelihood
   # (the sum over its rows) must match log_lik_binom() on the first category
diff --git a/tests/testthat/test-ran.R b/tests/testthat/test-ran.R
index 4b9886d6..1dd605e2 100644
--- a/tests/testthat/test-ran.R
+++ b/tests/testthat/test-ran.R
@@ -196,36 +196,25 @@ test_that("ran_multinom", {
     expect_identical(x[3] + x[4], 6L)
   })
 })
-test_that("ran_multinom with two categories matches ran_binom", {
-  size <- c(10, 4, 20)
-  prob <- c(0.2, 0.5, 0.75)
-  # rmultinom() draws a two-category trial with a single rbinom() call, so
-  # the first category is stream-identical to ran_binom()
-  withr::with_seed(42, {
-    x <- ran_multinom(
-      size = rep(size, each = 2),
-      prob = as.vector(rbind(prob, 1 - prob)),
-      group = rep(seq_along(size), each = 2)
-    )
-  })
-  withr::with_seed(42, {
-    y <- ran_binom(length(size), size = size, prob = prob)
-  })
-  expect_identical(matrix(x, nrow = 2)[1, ], y)
 
-  # and the marginal of a category is binomial, so its mean and variance
-  # match size * prob and size * prob * (1 - prob)
+test_that("ran_multinom with two categories draws a binomial first category", {
+  # a category count is marginally binomial, so over many trials the first
+  # category's mean, variance and count distribution match the binomial
+  size <- 10
+  prob <- 0.3
   n_group <- 20000
   withr::with_seed(7, {
     x <- ran_multinom(
-      size = rep(10, 2 * n_group),
-      prob = rep(c(0.3, 0.7), n_group),
+      size = rep(size, 2 * n_group),
+      prob = rep(c(prob, 1 - prob), n_group),
       group = rep(seq_len(n_group), each = 2)
     )
   })
   first <- matrix(x, nrow = 2)[1, ]
-  expect_equal(mean(first), 10 * 0.3, tolerance = 0.01)
-  expect_equal(var(first), 10 * 0.3 * 0.7, tolerance = 0.01)
+  expect_equal(mean(first), size * prob, tolerance = 0.01)
+  expect_equal(var(first), size * prob * (1 - prob), tolerance = 0.01)
+  proportion <- as.numeric(table(factor(first, levels = 0:size))) / n_group
+  expect_lt(max(abs(proportion - dbinom(0:size, size, prob))), 0.01)
 })
 
 test_that("ran_neg_binom", {
diff --git a/tests/testthat/test-res.R b/tests/testthat/test-res.R
index 4cc11645..1e3db93e 100644
--- a/tests/testthat/test-res.R
+++ b/tests/testthat/test-res.R
@@ -649,6 +649,7 @@ test_that("res_multinom simulate", {
     expect_equal(sd(res), 1.00164264126585)
   })
 })
+
 test_that("res_multinom with two categories matches res_binom", {
   x <- c(0, 3, 7, 10)
   size <- 10
@@ -677,23 +678,6 @@ test_that("res_multinom with two categories matches res_binom", {
     sign(res_type("dev")[1, ]),
     sign(res_binom(x, size, prob, type = "dev"))
   )
-
-  # a simulated two-category trial is a single rbinom() draw, so it is
-  # stream-identical to res_binom(simulate = TRUE)
-  withr::with_seed(3, {
-    sim <- res_multinom(
-      x_long,
-      size,
-      prob_long,
-      group,
-      type = "data",
-      simulate = TRUE
-    )
-  })
-  withr::with_seed(3, {
-    sim_binom <- res_binom(x, size, prob, type = "data", simulate = TRUE)
-  })
-  expect_equal(matrix(sim, nrow = 2)[1, ], sim_binom)
 })
 
 test_that("res_neg_binom", {

From a932e22f228037ce9d16a3f039f2b805c3f3a243 Mon Sep 17 00:00:00 2001
From: Joe Thorley <joe@poissonconsulting.ca>
Date: Tue, 25 Aug 2026 10:26:43 +0100
Subject: [PATCH 3/3] Replace the ran_multinom binomial check with a
 repeatability check

Assert that a seeded two-category call returns fixed values rather than
checking the marginal distribution over 20,000 trials.

Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
---
 tests/testthat/test-ran.R | 19 +++++--------------
 1 file changed, 5 insertions(+), 14 deletions(-)

diff --git a/tests/testthat/test-ran.R b/tests/testthat/test-ran.R
index 1dd605e2..ca1f98ea 100644
--- a/tests/testthat/test-ran.R
+++ b/tests/testthat/test-ran.R
@@ -197,24 +197,15 @@ test_that("ran_multinom", {
   })
 })
 
-test_that("ran_multinom with two categories draws a binomial first category", {
-  # a category count is marginally binomial, so over many trials the first
-  # category's mean, variance and count distribution match the binomial
-  size <- 10
-  prob <- 0.3
-  n_group <- 20000
+test_that("ran_multinom with two categories is repeatable", {
   withr::with_seed(7, {
     x <- ran_multinom(
-      size = rep(size, 2 * n_group),
-      prob = rep(c(prob, 1 - prob), n_group),
-      group = rep(seq_len(n_group), each = 2)
+      size = rep(10, 6),
+      prob = rep(c(0.3, 0.7), 3),
+      group = rep(1:3, each = 2)
     )
   })
-  first <- matrix(x, nrow = 2)[1, ]
-  expect_equal(mean(first), size * prob, tolerance = 0.01)
-  expect_equal(var(first), size * prob * (1 - prob), tolerance = 0.01)
-  proportion <- as.numeric(table(factor(first, levels = 0:size))) / n_group
-  expect_lt(max(abs(proportion - dbinom(0:size, size, prob))), 0.01)
+  expect_identical(x, c(6L, 4L, 3L, 7L, 1L, 9L))
 })
 
 test_that("ran_neg_binom", {