Test the multinomial functions against their binomial special case

tests/testthat/test-dev.R

@@ -570,6 +570,33 @@ test_that("dev_multinom res", {
   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(
     dev_neg_binom(integer(0), integer(0), integer(0)),

tests/testthat/test-log-lik.R

@@ -428,6 +428,24 @@ test_that("log_lik_multinom", {
   )
 })
 
+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(
     log_lik_neg_binom(0, 2, 1),

tests/testthat/test-ran.R

@@ -197,6 +197,17 @@ test_that("ran_multinom", {
   })
 })
 
+test_that("ran_multinom with two categories is repeatable", {
+  withr::with_seed(7, {
+    x <- ran_multinom(
+      size = rep(10, 6),
+      prob = rep(c(0.3, 0.7), 3),
+      group = rep(1:3, each = 2)
+    )
+  })
+  expect_identical(x, c(6L, 4L, 3L, 7L, 1L, 9L))
+})
+
 test_that("ran_neg_binom", {
   expect_error(ran_neg_binom(NA_integer_))
   expect_error(ran_neg_binom(integer(0)))

tests/testthat/test-res.R

@@ -650,6 +650,36 @@ test_that("res_multinom simulate", {
   })
 })
 
+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"))
+  )
+})
+
 test_that("res_neg_binom", {
   expect_identical(
     res_neg_binom(integer(0), integer(0), integer(0)),