From d33ad4a4a88bff3e74b062e2aaa4adf3d4d6f094 Mon Sep 17 00:00:00 2001 From: Nicole Hill Date: Wed, 5 Aug 2026 16:30:36 -0700 Subject: [PATCH 1/9] Add multinomial distribution functions (#62) Adds ran_multinom(), log_lik_multinom(), dev_multinom(), and res_multinom(). Data are expected in long format (one row per category per trial, identified by a `group` argument) so that residuals work one-per-observation with embr's posterior_predictive_check(), and so that a signed residual is always well defined per cell without needing an "extended sign" for a multi-category row. dev_multinom()/log_lik_multinom() rely on the identity that a multinomial is equivalent to independent Poissons conditional on the trial total, letting each category's contribution be computed from dev_pois()/log_lik_pois() alone. dev_multinom() is verified against glmnet's multinomial deviance() and log_lik_multinom() against dmultinom(). --- DESCRIPTION | 1 + NAMESPACE | 4 ++ R/dev.R | 34 ++++++++++++++++ R/internal.R | 22 ++++++++++ R/log-lik.R | 46 +++++++++++++++++++++ R/params.R | 4 ++ R/ran.R | 40 ++++++++++++++++++ R/res.R | 62 ++++++++++++++++++++++++++++ inst/WORDLIST | 1 + man/dev_bern.Rd | 1 + man/dev_beta_binom.Rd | 1 + man/dev_binom.Rd | 1 + man/dev_gamma.Rd | 1 + man/dev_gamma_pois.Rd | 1 + man/dev_lnorm.Rd | 1 + man/dev_multinom.Rd | 61 +++++++++++++++++++++++++++ man/dev_neg_binom.Rd | 1 + man/dev_norm.Rd | 1 + man/dev_pois.Rd | 1 + man/dev_pois_zi.Rd | 1 + man/dev_skewlnorm.Rd | 1 + man/dev_skewnorm.Rd | 1 + man/dev_student.Rd | 1 + man/log_lik_bern.Rd | 1 + man/log_lik_beta.Rd | 1 + man/log_lik_beta_binom.Rd | 1 + man/log_lik_binom.Rd | 1 + man/log_lik_exp.Rd | 1 + man/log_lik_gamma.Rd | 1 + man/log_lik_gamma_pois.Rd | 1 + man/log_lik_gamma_pois_zi.Rd | 1 + man/log_lik_lnorm.Rd | 1 + man/log_lik_multinom.Rd | 67 ++++++++++++++++++++++++++++++ man/log_lik_neg_binom.Rd | 1 + man/log_lik_norm.Rd | 1 + man/log_lik_pois.Rd | 1 + man/log_lik_pois_zi.Rd | 1 + man/log_lik_skewlnorm.Rd | 1 + man/log_lik_skewnorm.Rd | 1 + man/log_lik_student.Rd | 1 + man/log_lik_unif.Rd | 1 + man/params.Rd | 5 +++ man/ran_bern.Rd | 1 + man/ran_beta_binom.Rd | 1 + man/ran_binom.Rd | 1 + man/ran_gamma.Rd | 1 + man/ran_gamma_pois.Rd | 1 + man/ran_gamma_pois_zi.Rd | 1 + man/ran_lnorm.Rd | 1 + man/ran_multinom.Rd | 57 ++++++++++++++++++++++++++ man/ran_neg_binom.Rd | 1 + man/ran_norm.Rd | 1 + man/ran_pois.Rd | 1 + man/ran_pois_zi.Rd | 1 + man/ran_skewlnorm.Rd | 1 + man/ran_skewnorm.Rd | 1 + man/ran_student.Rd | 1 + man/res_bern.Rd | 1 + man/res_beta_binom.Rd | 1 + man/res_binom.Rd | 1 + man/res_gamma.Rd | 1 + man/res_gamma_pois.Rd | 1 + man/res_gamma_pois_zi.Rd | 1 + man/res_lnorm.Rd | 1 + man/res_multinom.Rd | 77 +++++++++++++++++++++++++++++++++++ man/res_neg_binom.Rd | 1 + man/res_norm.Rd | 1 + man/res_pois.Rd | 1 + man/res_pois_zi.Rd | 1 + man/res_skewlnorm.Rd | 1 + man/res_skewnorm.Rd | 1 + man/res_student.Rd | 1 + tests/testthat/test-dev.R | 73 +++++++++++++++++++++++++++++++++ tests/testthat/test-log-lik.R | 38 +++++++++++++++++ tests/testthat/test-ran.R | 29 +++++++++++++ tests/testthat/test-res.R | 39 ++++++++++++++++++ 76 files changed, 718 insertions(+) create mode 100644 man/dev_multinom.Rd create mode 100644 man/log_lik_multinom.Rd create mode 100644 man/ran_multinom.Rd create mode 100644 man/res_multinom.Rd diff --git a/DESCRIPTION b/DESCRIPTION index ad6aa191..3e7384ec 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -33,6 +33,7 @@ Suggests: covr, extraDistr, ggplot2, + glmnet, hms, knitr, memoise, diff --git a/NAMESPACE b/NAMESPACE index 1cb86856..0ac22bcf 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -39,6 +39,7 @@ export(dev_gamma) export(dev_gamma_pois) export(dev_gamma_pois_zi) export(dev_lnorm) +export(dev_multinom) export(dev_neg_binom) export(dev_norm) export(dev_pois) @@ -72,6 +73,7 @@ export(log_lik_gamma) export(log_lik_gamma_pois) export(log_lik_gamma_pois_zi) export(log_lik_lnorm) +export(log_lik_multinom) export(log_lik_neg_binom) export(log_lik_norm) export(log_lik_pois) @@ -150,6 +152,7 @@ export(ran_gamma) export(ran_gamma_pois) export(ran_gamma_pois_zi) export(ran_lnorm) +export(ran_multinom) export(ran_neg_binom) export(ran_norm) export(ran_pois) @@ -165,6 +168,7 @@ export(res_gamma) export(res_gamma_pois) export(res_gamma_pois_zi) export(res_lnorm) +export(res_multinom) export(res_neg_binom) export(res_norm) export(res_pois) diff --git a/R/dev.R b/R/dev.R index 641286c5..7bca89e7 100644 --- a/R/dev.R +++ b/R/dev.R @@ -267,6 +267,40 @@ dev_lnorm <- function(x, meanlog = 0, sdlog = 1, res = FALSE) { dev_norm(log(x), mean = meanlog, sd = sdlog, res = res) } +#' Multinomial Deviances +#' +#' The multinomial distribution models the counts across two or more +#' mutually exclusive categories arising from a fixed number of trials. Data +#' are in \emph{long} format: one row per category per trial (see +#' [log_lik_multinom()] for details of the `group` argument used elsewhere +#' in the multinomial family). +#' +#' Unlike the other `dev_*()` functions, `dev_multinom()` doesn't take a +#' `group` argument, because a category's deviance contribution only depends +#' on its own `x` and `mu = size * prob` -- the multinomial coefficient +#' cancels out of the deviance difference. This is the Poisson-equivalent +#' deviance (see [dev_pois()]): summing it over the rows of one trial +#' recovers the trial's exact multinomial deviance, provided `prob` sums to +#' 1 across those rows. Because each cell's deviance is an ordinary scalar +#' comparison of `x` to `mu`, its deviance residual has an ordinary sign -- +#' there's no need for the "extended sign" methods proposed elsewhere for +#' multinomial residuals. +#' +#' @inheritParams params +#' @param x A non-negative whole numeric vector of the category counts. +#' @param prob A numeric vector of the probability of the category. Must sum +#' to 1 across the rows belonging to the same trial. +#' +#' @return An numeric vector of the corresponding deviances or deviance residuals. +#' @family dev_dist +#' @export +#' +#' @examples +#' dev_multinom(c(1, 3, 6), size = 10, prob = c(0.2, 0.3, 0.5)) +dev_multinom <- function(x, size = 1, prob, res = FALSE) { + dev_pois(x, lambda = size * prob, res = res) +} + #' Negative Binomial Deviances #' #' @inheritParams params diff --git a/R/internal.R b/R/internal.R index a56e483b..9257b6a8 100644 --- a/R/internal.R +++ b/R/internal.R @@ -1,3 +1,25 @@ dev_res <- function(x, mu, dev) { sign(x - mu) * sqrt(dev) } + +# Validates that every row sharing a `group` (multinomial trial) has the +# same `size` and that its `prob` values sum to 1, since both the +# Poisson-equivalent deviance/log-likelihood identities and simulation via +# rmultinom() require this to hold. +chk_multinom_group <- function(size, prob, group) { + for (idx in split(seq_along(group), group)) { + if (length(unique(size[idx])) > 1L) { + stop( + "`size` must be the same for every row belonging to the same `group` (multinomial trial).", + call. = FALSE + ) + } + prob_sum <- sum(prob[idx]) + if (!is.na(prob_sum) && abs(prob_sum - 1) > 1e-6) { + stop( + "`prob` must sum to 1 for every `group` (multinomial trial).", + call. = FALSE + ) + } + } +} diff --git a/R/log-lik.R b/R/log-lik.R index c579fe64..eb83b436 100644 --- a/R/log-lik.R +++ b/R/log-lik.R @@ -427,6 +427,52 @@ log_lik_lnorm <- function(x, meanlog = 0, sdlog = 1, tlower = 0, tupper = Inf) { log_lik } +#' Multinomial Log-Likelihood +#' +#' The multinomial distribution models the counts across two or more +#' mutually exclusive categories arising from a fixed number of trials. Data +#' are in \emph{long} format: one row per category per trial, with `group` +#' identifying which rows belong to the same trial. All rows sharing a +#' `group` must have the same `size`, and their `prob` values must sum to 1. +#' +#' A trial's log-likelihood doesn't split evenly across its category rows, +#' because the multinomial coefficient is a property of the whole trial, not +#' any one category. `log_lik_multinom()` uses the identity that a +#' multinomial is equivalent to independent Poissons conditional on the +#' trial total: the log-likelihood of category `k` of trial `i` is the +#' Poisson log-likelihood of `x` given `mu = size * prob`, minus an even +#' share of the trial's normalizing constant (so that summing +#' `log_lik_multinom()` over the rows of one `group` recovers the trial's +#' exact multinomial log-likelihood). +#' +#' @inheritParams params +#' @param x A non-negative whole numeric vector of the category counts. +#' @param prob A numeric vector of the probability of the category. Must sum +#' to 1 across the rows sharing the same `group`. +#' +#' @return An numeric vector of the corresponding log-likelihoods, one value +#' per row of `x`. +#' @family log_lik_dist +#' @export +#' +#' @examples +#' log_lik_multinom(c(1, 3, 6), size = 10, prob = c(0.2, 0.3, 0.5), group = c(1, 1, 1)) +log_lik_multinom <- function(x, size = 1, prob, group) { + n <- length(x) + if (!n || !length(size) || !length(prob) || !length(group)) { + return(numeric(0)) + } + size <- rep_len(size, n) + prob <- rep_len(prob, n) + group <- rep_len(group, n) + chk_multinom_group(size, prob, group) + mu <- size * prob + log_lik <- log_lik_pois(x, mu) + k <- ave(seq_along(group), group, FUN = length) + const <- log_lik_pois(size, size) + log_lik - const / k +} + #' Negative Binomial Log-Likelihood #' #' @inheritParams params diff --git a/R/params.R b/R/params.R index ec79aed9..513fb8a5 100644 --- a/R/params.R +++ b/R/params.R @@ -16,6 +16,10 @@ #' level. #' @param directional A flag specifying whether probabilities less than 0.5 #' should be returned as negative values. +#' @param group A vector identifying which rows belong to the same +#' multinomial trial, i.e., whose `x` values must sum to `size` and whose +#' `prob` values must sum to 1. Required because the categories within a +#' single multinomial trial are not independent. #' @param lambda A non-negative numeric vector of means. #' @param level A number > 0 and <= 1 specifying the probability coverage of the #' interval. diff --git a/R/ran.R b/R/ran.R index 65369c77..dbd963c5 100644 --- a/R/ran.R +++ b/R/ran.R @@ -113,6 +113,46 @@ ran_lnorm <- function(n = 1, meanlog = 0, sdlog = 1) { stats::rlnorm(n, meanlog = meanlog, sdlog = sdlog) } +#' Multinomial Random Samples +#' +#' The multinomial distribution models the counts across two or more +#' mutually exclusive categories arising from a fixed number of trials. Data +#' (and therefore random samples) are in \emph{long} format: one value per +#' category per trial, with `group` identifying which rows belong to the +#' same trial. All rows sharing a `group` must have the same `size`, and +#' their `prob` values must sum to 1. +#' +#' Unlike the other `ran_*()` functions, `ran_multinom()` has no `n` +#' argument: the number of samples is fully determined by `length(prob)` +#' (equivalently `length(group)`), because a trial's categories can't be +#' generated independently of one another. +#' +#' @inheritParams params +#' @param prob A numeric vector of the probability of the category. Must sum +#' to 1 across the rows sharing the same `group`. +#' @return An integer vector of the random samples, one per row of `prob`. +#' @family ran_dist +#' @export +#' +#' @examples +#' ran_multinom(size = 10, prob = c(0.2, 0.3, 0.5), group = c(1, 1, 1)) +ran_multinom <- function(size = 1, prob, group) { + n <- length(prob) + if (!n) { + return(integer(0)) + } + chk_compatible_lengths(rep(1, n), size, group) + size <- rep_len(size, n) + prob <- rep_len(prob, n) + group <- rep_len(group, n) + chk_multinom_group(size, prob, group) + x <- rep(NA_real_, n) + for (idx in split(seq_len(n), group)) { + x[idx] <- stats::rmultinom(1, size = size[idx[1]], prob = prob[idx])[, 1] + } + as.integer(x) +} + #' Negative Binomial Random Samples #' #' Identical to Gamma-Poisson Random Samples. diff --git a/R/res.R b/R/res.R index c260d5e5..2b701063 100644 --- a/R/res.R +++ b/R/res.R @@ -251,6 +251,68 @@ res_lnorm <- function( ) } +#' Multinomial Residuals +#' +#' The multinomial distribution models the counts across two or more +#' mutually exclusive categories arising from a fixed number of trials. Data +#' are in \emph{long} format: one row per category per trial, with `group` +#' identifying which rows belong to the same trial (see +#' [log_lik_multinom()] for details). `res_multinom()` returns one residual +#' per row (per category per trial), not one per trial -- since a trial's +#' categories are not independent, there's no single meaningful residual (or +#' sign) for the trial as a whole. The classic per-trial deviance statistic +#' can always be recovered by summing the squared `type = "dev"` residuals +#' within a `group`. +#' +#' Because a category count is marginally binomial (`x ~ Binomial(size, +#' prob)`), the `"standardized"` residual uses the same formula as +#' [res_binom()]; this is already the covariance-adjusted (unit-variance) +#' residual, not the naive independent-Poisson version. +#' +#' `group` is only used when `simulate = TRUE`, to draw a joint, +#' correlation-preserving replicate for each trial (via [ran_multinom()]) +#' rather than simulating each category independently, which would +#' understate the true (negative) covariance among a trial's categories. For +#' the simulated residuals from this function to be valid inputs to +#' `embr::posterior_predictive_check()`, the model must be fit with +#' `new_expr_vec = TRUE`, so that `res_multinom()` is evaluated once on the +#' full data vector (with visibility of every row in each `group`) rather +#' than row by row. +#' +#' @inheritParams params +#' @param x A non-negative whole numeric vector of the category counts. +#' @param prob A numeric vector of the probability of the category. Must sum +#' to 1 across the rows sharing the same `group`. +#' +#' @return An numeric vector of the corresponding residuals. +#' @family res_dist +#' @export +#' +#' @examples +#' res_multinom(c(1, 3, 6), size = 10, prob = c(0.2, 0.3, 0.5), group = c(1, 1, 1)) +res_multinom <- function( + x, + size = 1, + prob, + group, + type = "dev", + simulate = FALSE +) { + chk_string(type) + if (!vld_false(simulate)) { + x <- ran_multinom(size = size, prob = prob, group = group) + } + mu <- size * prob + switch( + type, + data = x, + raw = x - mu, + standardized = (x - mu) / sqrt(mu * (1 - prob)), + dev = dev_multinom(x, size = size, prob = prob, res = TRUE), + chk_subset(x, c("data", "raw", "dev", "standardized")) + ) +} + #' Negative Binomial Residuals #' #' @inheritParams params diff --git a/inst/WORDLIST b/inst/WORDLIST index 127f2b03..8257c2a9 100644 --- a/inst/WORDLIST +++ b/inst/WORDLIST @@ -11,6 +11,7 @@ Numericise Numericize ORCID POSIXct +Poissons Psychol Schaub Shachar diff --git a/man/dev_bern.Rd b/man/dev_bern.Rd index f978015f..91e829b8 100644 --- a/man/dev_bern.Rd +++ b/man/dev_bern.Rd @@ -31,6 +31,7 @@ Other dev_dist: \code{\link[=dev_gamma]{dev_gamma()}}, \code{\link[=dev_gamma_pois]{dev_gamma_pois()}}, \code{\link[=dev_lnorm]{dev_lnorm()}}, +\code{\link[=dev_multinom]{dev_multinom()}}, \code{\link[=dev_neg_binom]{dev_neg_binom()}}, \code{\link[=dev_norm]{dev_norm()}}, \code{\link[=dev_pois]{dev_pois()}}, diff --git a/man/dev_beta_binom.Rd b/man/dev_beta_binom.Rd index 84ca558c..7edf8110 100644 --- a/man/dev_beta_binom.Rd +++ b/man/dev_beta_binom.Rd @@ -45,6 +45,7 @@ Other dev_dist: \code{\link[=dev_gamma]{dev_gamma()}}, \code{\link[=dev_gamma_pois]{dev_gamma_pois()}}, \code{\link[=dev_lnorm]{dev_lnorm()}}, +\code{\link[=dev_multinom]{dev_multinom()}}, \code{\link[=dev_neg_binom]{dev_neg_binom()}}, \code{\link[=dev_norm]{dev_norm()}}, \code{\link[=dev_pois]{dev_pois()}}, diff --git a/man/dev_binom.Rd b/man/dev_binom.Rd index 4d83f953..063a09fb 100644 --- a/man/dev_binom.Rd +++ b/man/dev_binom.Rd @@ -33,6 +33,7 @@ Other dev_dist: \code{\link[=dev_gamma]{dev_gamma()}}, \code{\link[=dev_gamma_pois]{dev_gamma_pois()}}, \code{\link[=dev_lnorm]{dev_lnorm()}}, +\code{\link[=dev_multinom]{dev_multinom()}}, \code{\link[=dev_neg_binom]{dev_neg_binom()}}, \code{\link[=dev_norm]{dev_norm()}}, \code{\link[=dev_pois]{dev_pois()}}, diff --git a/man/dev_gamma.Rd b/man/dev_gamma.Rd index 8041ccfe..cf9ebf4d 100644 --- a/man/dev_gamma.Rd +++ b/man/dev_gamma.Rd @@ -32,6 +32,7 @@ Other dev_dist: \code{\link[=dev_binom]{dev_binom()}}, \code{\link[=dev_gamma_pois]{dev_gamma_pois()}}, \code{\link[=dev_lnorm]{dev_lnorm()}}, +\code{\link[=dev_multinom]{dev_multinom()}}, \code{\link[=dev_neg_binom]{dev_neg_binom()}}, \code{\link[=dev_norm]{dev_norm()}}, \code{\link[=dev_pois]{dev_pois()}}, diff --git a/man/dev_gamma_pois.Rd b/man/dev_gamma_pois.Rd index 89f97a4f..1f68251b 100644 --- a/man/dev_gamma_pois.Rd +++ b/man/dev_gamma_pois.Rd @@ -33,6 +33,7 @@ Other dev_dist: \code{\link[=dev_binom]{dev_binom()}}, \code{\link[=dev_gamma]{dev_gamma()}}, \code{\link[=dev_lnorm]{dev_lnorm()}}, +\code{\link[=dev_multinom]{dev_multinom()}}, \code{\link[=dev_neg_binom]{dev_neg_binom()}}, \code{\link[=dev_norm]{dev_norm()}}, \code{\link[=dev_pois]{dev_pois()}}, diff --git a/man/dev_lnorm.Rd b/man/dev_lnorm.Rd index 0f10c26c..6512b9d3 100644 --- a/man/dev_lnorm.Rd +++ b/man/dev_lnorm.Rd @@ -33,6 +33,7 @@ Other dev_dist: \code{\link[=dev_binom]{dev_binom()}}, \code{\link[=dev_gamma]{dev_gamma()}}, \code{\link[=dev_gamma_pois]{dev_gamma_pois()}}, +\code{\link[=dev_multinom]{dev_multinom()}}, \code{\link[=dev_neg_binom]{dev_neg_binom()}}, \code{\link[=dev_norm]{dev_norm()}}, \code{\link[=dev_pois]{dev_pois()}}, diff --git a/man/dev_multinom.Rd b/man/dev_multinom.Rd new file mode 100644 index 00000000..472ac5e1 --- /dev/null +++ b/man/dev_multinom.Rd @@ -0,0 +1,61 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/dev.R +\name{dev_multinom} +\alias{dev_multinom} +\title{Multinomial Deviances} +\usage{ +dev_multinom(x, size = 1, prob, res = FALSE) +} +\arguments{ +\item{x}{A non-negative whole numeric vector of the category counts.} + +\item{size}{A non-negative whole numeric vector of the number of trials.} + +\item{prob}{A numeric vector of the probability of the category. Must sum +to 1 across the rows belonging to the same trial.} + +\item{res}{A flag specifying whether to return the deviance residual as +opposed to the deviance.} +} +\value{ +An numeric vector of the corresponding deviances or deviance residuals. +} +\description{ +The multinomial distribution models the counts across two or more +mutually exclusive categories arising from a fixed number of trials. Data +are in \emph{long} format: one row per category per trial (see +\code{\link[=log_lik_multinom]{log_lik_multinom()}} for details of the \code{group} argument used elsewhere +in the multinomial family). +} +\details{ +Unlike the other \verb{dev_*()} functions, \code{dev_multinom()} doesn't take a +\code{group} argument, because a category's deviance contribution only depends +on its own \code{x} and \code{mu = size * prob} -- the multinomial coefficient +cancels out of the deviance difference. This is the Poisson-equivalent +deviance (see \code{\link[=dev_pois]{dev_pois()}}): summing it over the rows of one trial +recovers the trial's exact multinomial deviance, provided \code{prob} sums to +1 across those rows. Because each cell's deviance is an ordinary scalar +comparison of \code{x} to \code{mu}, its deviance residual has an ordinary sign -- +there's no need for the "extended sign" methods proposed elsewhere for +multinomial residuals. +} +\examples{ +dev_multinom(c(1, 3, 6), size = 10, prob = c(0.2, 0.3, 0.5)) +} +\seealso{ +Other dev_dist: +\code{\link[=dev_bern]{dev_bern()}}, +\code{\link[=dev_beta_binom]{dev_beta_binom()}}, +\code{\link[=dev_binom]{dev_binom()}}, +\code{\link[=dev_gamma]{dev_gamma()}}, +\code{\link[=dev_gamma_pois]{dev_gamma_pois()}}, +\code{\link[=dev_lnorm]{dev_lnorm()}}, +\code{\link[=dev_neg_binom]{dev_neg_binom()}}, +\code{\link[=dev_norm]{dev_norm()}}, +\code{\link[=dev_pois]{dev_pois()}}, +\code{\link[=dev_pois_zi]{dev_pois_zi()}}, +\code{\link[=dev_skewlnorm]{dev_skewlnorm()}}, +\code{\link[=dev_skewnorm]{dev_skewnorm()}}, +\code{\link[=dev_student]{dev_student()}} +} +\concept{dev_dist} diff --git a/man/dev_neg_binom.Rd b/man/dev_neg_binom.Rd index 5696099c..947dcf15 100644 --- a/man/dev_neg_binom.Rd +++ b/man/dev_neg_binom.Rd @@ -34,6 +34,7 @@ Other dev_dist: \code{\link[=dev_gamma]{dev_gamma()}}, \code{\link[=dev_gamma_pois]{dev_gamma_pois()}}, \code{\link[=dev_lnorm]{dev_lnorm()}}, +\code{\link[=dev_multinom]{dev_multinom()}}, \code{\link[=dev_norm]{dev_norm()}}, \code{\link[=dev_pois]{dev_pois()}}, \code{\link[=dev_pois_zi]{dev_pois_zi()}}, diff --git a/man/dev_norm.Rd b/man/dev_norm.Rd index 5b556302..fc685b14 100644 --- a/man/dev_norm.Rd +++ b/man/dev_norm.Rd @@ -33,6 +33,7 @@ Other dev_dist: \code{\link[=dev_gamma]{dev_gamma()}}, \code{\link[=dev_gamma_pois]{dev_gamma_pois()}}, \code{\link[=dev_lnorm]{dev_lnorm()}}, +\code{\link[=dev_multinom]{dev_multinom()}}, \code{\link[=dev_neg_binom]{dev_neg_binom()}}, \code{\link[=dev_pois]{dev_pois()}}, \code{\link[=dev_pois_zi]{dev_pois_zi()}}, diff --git a/man/dev_pois.Rd b/man/dev_pois.Rd index b99db9cd..f00cd65a 100644 --- a/man/dev_pois.Rd +++ b/man/dev_pois.Rd @@ -31,6 +31,7 @@ Other dev_dist: \code{\link[=dev_gamma]{dev_gamma()}}, \code{\link[=dev_gamma_pois]{dev_gamma_pois()}}, \code{\link[=dev_lnorm]{dev_lnorm()}}, +\code{\link[=dev_multinom]{dev_multinom()}}, \code{\link[=dev_neg_binom]{dev_neg_binom()}}, \code{\link[=dev_norm]{dev_norm()}}, \code{\link[=dev_pois_zi]{dev_pois_zi()}}, diff --git a/man/dev_pois_zi.Rd b/man/dev_pois_zi.Rd index 23db2ccf..7aa37c09 100644 --- a/man/dev_pois_zi.Rd +++ b/man/dev_pois_zi.Rd @@ -34,6 +34,7 @@ Other dev_dist: \code{\link[=dev_gamma]{dev_gamma()}}, \code{\link[=dev_gamma_pois]{dev_gamma_pois()}}, \code{\link[=dev_lnorm]{dev_lnorm()}}, +\code{\link[=dev_multinom]{dev_multinom()}}, \code{\link[=dev_neg_binom]{dev_neg_binom()}}, \code{\link[=dev_norm]{dev_norm()}}, \code{\link[=dev_pois]{dev_pois()}}, diff --git a/man/dev_skewlnorm.Rd b/man/dev_skewlnorm.Rd index 5b39cef8..c3148042 100644 --- a/man/dev_skewlnorm.Rd +++ b/man/dev_skewlnorm.Rd @@ -40,6 +40,7 @@ Other dev_dist: \code{\link[=dev_gamma]{dev_gamma()}}, \code{\link[=dev_gamma_pois]{dev_gamma_pois()}}, \code{\link[=dev_lnorm]{dev_lnorm()}}, +\code{\link[=dev_multinom]{dev_multinom()}}, \code{\link[=dev_neg_binom]{dev_neg_binom()}}, \code{\link[=dev_norm]{dev_norm()}}, \code{\link[=dev_pois]{dev_pois()}}, diff --git a/man/dev_skewnorm.Rd b/man/dev_skewnorm.Rd index 389cea5d..c46bd9b7 100644 --- a/man/dev_skewnorm.Rd +++ b/man/dev_skewnorm.Rd @@ -39,6 +39,7 @@ Other dev_dist: \code{\link[=dev_gamma]{dev_gamma()}}, \code{\link[=dev_gamma_pois]{dev_gamma_pois()}}, \code{\link[=dev_lnorm]{dev_lnorm()}}, +\code{\link[=dev_multinom]{dev_multinom()}}, \code{\link[=dev_neg_binom]{dev_neg_binom()}}, \code{\link[=dev_norm]{dev_norm()}}, \code{\link[=dev_pois]{dev_pois()}}, diff --git a/man/dev_student.Rd b/man/dev_student.Rd index 3f9739d6..6ce0ec43 100644 --- a/man/dev_student.Rd +++ b/man/dev_student.Rd @@ -36,6 +36,7 @@ Other dev_dist: \code{\link[=dev_gamma]{dev_gamma()}}, \code{\link[=dev_gamma_pois]{dev_gamma_pois()}}, \code{\link[=dev_lnorm]{dev_lnorm()}}, +\code{\link[=dev_multinom]{dev_multinom()}}, \code{\link[=dev_neg_binom]{dev_neg_binom()}}, \code{\link[=dev_norm]{dev_norm()}}, \code{\link[=dev_pois]{dev_pois()}}, diff --git a/man/log_lik_bern.Rd b/man/log_lik_bern.Rd index 18ddb194..9212e2ae 100644 --- a/man/log_lik_bern.Rd +++ b/man/log_lik_bern.Rd @@ -31,6 +31,7 @@ Other log_lik_dist: \code{\link[=log_lik_gamma_pois]{log_lik_gamma_pois()}}, \code{\link[=log_lik_gamma_pois_zi]{log_lik_gamma_pois_zi()}}, \code{\link[=log_lik_lnorm]{log_lik_lnorm()}}, +\code{\link[=log_lik_multinom]{log_lik_multinom()}}, \code{\link[=log_lik_neg_binom]{log_lik_neg_binom()}}, \code{\link[=log_lik_norm]{log_lik_norm()}}, \code{\link[=log_lik_pois]{log_lik_pois()}}, diff --git a/man/log_lik_beta.Rd b/man/log_lik_beta.Rd index 9e6453fb..d4def7b0 100644 --- a/man/log_lik_beta.Rd +++ b/man/log_lik_beta.Rd @@ -36,6 +36,7 @@ Other log_lik_dist: \code{\link[=log_lik_gamma_pois]{log_lik_gamma_pois()}}, \code{\link[=log_lik_gamma_pois_zi]{log_lik_gamma_pois_zi()}}, \code{\link[=log_lik_lnorm]{log_lik_lnorm()}}, +\code{\link[=log_lik_multinom]{log_lik_multinom()}}, \code{\link[=log_lik_neg_binom]{log_lik_neg_binom()}}, \code{\link[=log_lik_norm]{log_lik_norm()}}, \code{\link[=log_lik_pois]{log_lik_pois()}}, diff --git a/man/log_lik_beta_binom.Rd b/man/log_lik_beta_binom.Rd index 8c737976..8528b741 100644 --- a/man/log_lik_beta_binom.Rd +++ b/man/log_lik_beta_binom.Rd @@ -58,6 +58,7 @@ Other log_lik_dist: \code{\link[=log_lik_gamma_pois]{log_lik_gamma_pois()}}, \code{\link[=log_lik_gamma_pois_zi]{log_lik_gamma_pois_zi()}}, \code{\link[=log_lik_lnorm]{log_lik_lnorm()}}, +\code{\link[=log_lik_multinom]{log_lik_multinom()}}, \code{\link[=log_lik_neg_binom]{log_lik_neg_binom()}}, \code{\link[=log_lik_norm]{log_lik_norm()}}, \code{\link[=log_lik_pois]{log_lik_pois()}}, diff --git a/man/log_lik_binom.Rd b/man/log_lik_binom.Rd index 475d74f8..5e90cf80 100644 --- a/man/log_lik_binom.Rd +++ b/man/log_lik_binom.Rd @@ -37,6 +37,7 @@ Other log_lik_dist: \code{\link[=log_lik_gamma_pois]{log_lik_gamma_pois()}}, \code{\link[=log_lik_gamma_pois_zi]{log_lik_gamma_pois_zi()}}, \code{\link[=log_lik_lnorm]{log_lik_lnorm()}}, +\code{\link[=log_lik_multinom]{log_lik_multinom()}}, \code{\link[=log_lik_neg_binom]{log_lik_neg_binom()}}, \code{\link[=log_lik_norm]{log_lik_norm()}}, \code{\link[=log_lik_pois]{log_lik_pois()}}, diff --git a/man/log_lik_exp.Rd b/man/log_lik_exp.Rd index 5192ae8f..94f16a05 100644 --- a/man/log_lik_exp.Rd +++ b/man/log_lik_exp.Rd @@ -34,6 +34,7 @@ Other log_lik_dist: \code{\link[=log_lik_gamma_pois]{log_lik_gamma_pois()}}, \code{\link[=log_lik_gamma_pois_zi]{log_lik_gamma_pois_zi()}}, \code{\link[=log_lik_lnorm]{log_lik_lnorm()}}, +\code{\link[=log_lik_multinom]{log_lik_multinom()}}, \code{\link[=log_lik_neg_binom]{log_lik_neg_binom()}}, \code{\link[=log_lik_norm]{log_lik_norm()}}, \code{\link[=log_lik_pois]{log_lik_pois()}}, diff --git a/man/log_lik_gamma.Rd b/man/log_lik_gamma.Rd index 953dfd9e..9148ff42 100644 --- a/man/log_lik_gamma.Rd +++ b/man/log_lik_gamma.Rd @@ -36,6 +36,7 @@ Other log_lik_dist: \code{\link[=log_lik_gamma_pois]{log_lik_gamma_pois()}}, \code{\link[=log_lik_gamma_pois_zi]{log_lik_gamma_pois_zi()}}, \code{\link[=log_lik_lnorm]{log_lik_lnorm()}}, +\code{\link[=log_lik_multinom]{log_lik_multinom()}}, \code{\link[=log_lik_neg_binom]{log_lik_neg_binom()}}, \code{\link[=log_lik_norm]{log_lik_norm()}}, \code{\link[=log_lik_pois]{log_lik_pois()}}, diff --git a/man/log_lik_gamma_pois.Rd b/man/log_lik_gamma_pois.Rd index 51f94e71..47e1c63b 100644 --- a/man/log_lik_gamma_pois.Rd +++ b/man/log_lik_gamma_pois.Rd @@ -37,6 +37,7 @@ Other log_lik_dist: \code{\link[=log_lik_gamma]{log_lik_gamma()}}, \code{\link[=log_lik_gamma_pois_zi]{log_lik_gamma_pois_zi()}}, \code{\link[=log_lik_lnorm]{log_lik_lnorm()}}, +\code{\link[=log_lik_multinom]{log_lik_multinom()}}, \code{\link[=log_lik_neg_binom]{log_lik_neg_binom()}}, \code{\link[=log_lik_norm]{log_lik_norm()}}, \code{\link[=log_lik_pois]{log_lik_pois()}}, diff --git a/man/log_lik_gamma_pois_zi.Rd b/man/log_lik_gamma_pois_zi.Rd index 1700d074..1f8e102d 100644 --- a/man/log_lik_gamma_pois_zi.Rd +++ b/man/log_lik_gamma_pois_zi.Rd @@ -47,6 +47,7 @@ Other log_lik_dist: \code{\link[=log_lik_gamma]{log_lik_gamma()}}, \code{\link[=log_lik_gamma_pois]{log_lik_gamma_pois()}}, \code{\link[=log_lik_lnorm]{log_lik_lnorm()}}, +\code{\link[=log_lik_multinom]{log_lik_multinom()}}, \code{\link[=log_lik_neg_binom]{log_lik_neg_binom()}}, \code{\link[=log_lik_norm]{log_lik_norm()}}, \code{\link[=log_lik_pois]{log_lik_pois()}}, diff --git a/man/log_lik_lnorm.Rd b/man/log_lik_lnorm.Rd index d02d6233..ac7e102f 100644 --- a/man/log_lik_lnorm.Rd +++ b/man/log_lik_lnorm.Rd @@ -37,6 +37,7 @@ Other log_lik_dist: \code{\link[=log_lik_gamma]{log_lik_gamma()}}, \code{\link[=log_lik_gamma_pois]{log_lik_gamma_pois()}}, \code{\link[=log_lik_gamma_pois_zi]{log_lik_gamma_pois_zi()}}, +\code{\link[=log_lik_multinom]{log_lik_multinom()}}, \code{\link[=log_lik_neg_binom]{log_lik_neg_binom()}}, \code{\link[=log_lik_norm]{log_lik_norm()}}, \code{\link[=log_lik_pois]{log_lik_pois()}}, diff --git a/man/log_lik_multinom.Rd b/man/log_lik_multinom.Rd new file mode 100644 index 00000000..65fc845d --- /dev/null +++ b/man/log_lik_multinom.Rd @@ -0,0 +1,67 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/log-lik.R +\name{log_lik_multinom} +\alias{log_lik_multinom} +\title{Multinomial Log-Likelihood} +\usage{ +log_lik_multinom(x, size = 1, prob, group) +} +\arguments{ +\item{x}{A non-negative whole numeric vector of the category counts.} + +\item{size}{A non-negative whole numeric vector of the number of trials.} + +\item{prob}{A numeric vector of the probability of the category. Must sum +to 1 across the rows sharing the same \code{group}.} + +\item{group}{A vector identifying which rows belong to the same +multinomial trial, i.e., whose \code{x} values must sum to \code{size} and whose +\code{prob} values must sum to 1. Required because the categories within a +single multinomial trial are not independent.} +} +\value{ +An numeric vector of the corresponding log-likelihoods, one value +per row of \code{x}. +} +\description{ +The multinomial distribution models the counts across two or more +mutually exclusive categories arising from a fixed number of trials. Data +are in \emph{long} format: one row per category per trial, with \code{group} +identifying which rows belong to the same trial. All rows sharing a +\code{group} must have the same \code{size}, and their \code{prob} values must sum to 1. +} +\details{ +A trial's log-likelihood doesn't split evenly across its category rows, +because the multinomial coefficient is a property of the whole trial, not +any one category. \code{log_lik_multinom()} uses the identity that a +multinomial is equivalent to independent Poissons conditional on the +trial total: the log-likelihood of category \code{k} of trial \code{i} is the +Poisson log-likelihood of \code{x} given \code{mu = size * prob}, minus an even +share of the trial's normalizing constant (so that summing +\code{log_lik_multinom()} over the rows of one \code{group} recovers the trial's +exact multinomial log-likelihood). +} +\examples{ +log_lik_multinom(c(1, 3, 6), size = 10, prob = c(0.2, 0.3, 0.5), group = c(1, 1, 1)) +} +\seealso{ +Other log_lik_dist: +\code{\link[=log_lik_bern]{log_lik_bern()}}, +\code{\link[=log_lik_beta]{log_lik_beta()}}, +\code{\link[=log_lik_beta_binom]{log_lik_beta_binom()}}, +\code{\link[=log_lik_binom]{log_lik_binom()}}, +\code{\link[=log_lik_exp]{log_lik_exp()}}, +\code{\link[=log_lik_gamma]{log_lik_gamma()}}, +\code{\link[=log_lik_gamma_pois]{log_lik_gamma_pois()}}, +\code{\link[=log_lik_gamma_pois_zi]{log_lik_gamma_pois_zi()}}, +\code{\link[=log_lik_lnorm]{log_lik_lnorm()}}, +\code{\link[=log_lik_neg_binom]{log_lik_neg_binom()}}, +\code{\link[=log_lik_norm]{log_lik_norm()}}, +\code{\link[=log_lik_pois]{log_lik_pois()}}, +\code{\link[=log_lik_pois_zi]{log_lik_pois_zi()}}, +\code{\link[=log_lik_skewlnorm]{log_lik_skewlnorm()}}, +\code{\link[=log_lik_skewnorm]{log_lik_skewnorm()}}, +\code{\link[=log_lik_student]{log_lik_student()}}, +\code{\link[=log_lik_unif]{log_lik_unif()}} +} +\concept{log_lik_dist} diff --git a/man/log_lik_neg_binom.Rd b/man/log_lik_neg_binom.Rd index a4e7609c..c06b5259 100644 --- a/man/log_lik_neg_binom.Rd +++ b/man/log_lik_neg_binom.Rd @@ -38,6 +38,7 @@ Other log_lik_dist: \code{\link[=log_lik_gamma_pois]{log_lik_gamma_pois()}}, \code{\link[=log_lik_gamma_pois_zi]{log_lik_gamma_pois_zi()}}, \code{\link[=log_lik_lnorm]{log_lik_lnorm()}}, +\code{\link[=log_lik_multinom]{log_lik_multinom()}}, \code{\link[=log_lik_norm]{log_lik_norm()}}, \code{\link[=log_lik_pois]{log_lik_pois()}}, \code{\link[=log_lik_pois_zi]{log_lik_pois_zi()}}, diff --git a/man/log_lik_norm.Rd b/man/log_lik_norm.Rd index 8018e102..feed9c0a 100644 --- a/man/log_lik_norm.Rd +++ b/man/log_lik_norm.Rd @@ -37,6 +37,7 @@ Other log_lik_dist: \code{\link[=log_lik_gamma_pois]{log_lik_gamma_pois()}}, \code{\link[=log_lik_gamma_pois_zi]{log_lik_gamma_pois_zi()}}, \code{\link[=log_lik_lnorm]{log_lik_lnorm()}}, +\code{\link[=log_lik_multinom]{log_lik_multinom()}}, \code{\link[=log_lik_neg_binom]{log_lik_neg_binom()}}, \code{\link[=log_lik_pois]{log_lik_pois()}}, \code{\link[=log_lik_pois_zi]{log_lik_pois_zi()}}, diff --git a/man/log_lik_pois.Rd b/man/log_lik_pois.Rd index 0eb155e4..421c67ff 100644 --- a/man/log_lik_pois.Rd +++ b/man/log_lik_pois.Rd @@ -35,6 +35,7 @@ Other log_lik_dist: \code{\link[=log_lik_gamma_pois]{log_lik_gamma_pois()}}, \code{\link[=log_lik_gamma_pois_zi]{log_lik_gamma_pois_zi()}}, \code{\link[=log_lik_lnorm]{log_lik_lnorm()}}, +\code{\link[=log_lik_multinom]{log_lik_multinom()}}, \code{\link[=log_lik_neg_binom]{log_lik_neg_binom()}}, \code{\link[=log_lik_norm]{log_lik_norm()}}, \code{\link[=log_lik_pois_zi]{log_lik_pois_zi()}}, diff --git a/man/log_lik_pois_zi.Rd b/man/log_lik_pois_zi.Rd index b1e67d4b..8138a82f 100644 --- a/man/log_lik_pois_zi.Rd +++ b/man/log_lik_pois_zi.Rd @@ -38,6 +38,7 @@ Other log_lik_dist: \code{\link[=log_lik_gamma_pois]{log_lik_gamma_pois()}}, \code{\link[=log_lik_gamma_pois_zi]{log_lik_gamma_pois_zi()}}, \code{\link[=log_lik_lnorm]{log_lik_lnorm()}}, +\code{\link[=log_lik_multinom]{log_lik_multinom()}}, \code{\link[=log_lik_neg_binom]{log_lik_neg_binom()}}, \code{\link[=log_lik_norm]{log_lik_norm()}}, \code{\link[=log_lik_pois]{log_lik_pois()}}, diff --git a/man/log_lik_skewlnorm.Rd b/man/log_lik_skewlnorm.Rd index 99660b49..d9a9fbbb 100644 --- a/man/log_lik_skewlnorm.Rd +++ b/man/log_lik_skewlnorm.Rd @@ -51,6 +51,7 @@ Other log_lik_dist: \code{\link[=log_lik_gamma_pois]{log_lik_gamma_pois()}}, \code{\link[=log_lik_gamma_pois_zi]{log_lik_gamma_pois_zi()}}, \code{\link[=log_lik_lnorm]{log_lik_lnorm()}}, +\code{\link[=log_lik_multinom]{log_lik_multinom()}}, \code{\link[=log_lik_neg_binom]{log_lik_neg_binom()}}, \code{\link[=log_lik_norm]{log_lik_norm()}}, \code{\link[=log_lik_pois]{log_lik_pois()}}, diff --git a/man/log_lik_skewnorm.Rd b/man/log_lik_skewnorm.Rd index f4acf7f7..ce0a9f2a 100644 --- a/man/log_lik_skewnorm.Rd +++ b/man/log_lik_skewnorm.Rd @@ -43,6 +43,7 @@ Other log_lik_dist: \code{\link[=log_lik_gamma_pois]{log_lik_gamma_pois()}}, \code{\link[=log_lik_gamma_pois_zi]{log_lik_gamma_pois_zi()}}, \code{\link[=log_lik_lnorm]{log_lik_lnorm()}}, +\code{\link[=log_lik_multinom]{log_lik_multinom()}}, \code{\link[=log_lik_neg_binom]{log_lik_neg_binom()}}, \code{\link[=log_lik_norm]{log_lik_norm()}}, \code{\link[=log_lik_pois]{log_lik_pois()}}, diff --git a/man/log_lik_student.Rd b/man/log_lik_student.Rd index 31e83fa1..d3700d73 100644 --- a/man/log_lik_student.Rd +++ b/man/log_lik_student.Rd @@ -40,6 +40,7 @@ Other log_lik_dist: \code{\link[=log_lik_gamma_pois]{log_lik_gamma_pois()}}, \code{\link[=log_lik_gamma_pois_zi]{log_lik_gamma_pois_zi()}}, \code{\link[=log_lik_lnorm]{log_lik_lnorm()}}, +\code{\link[=log_lik_multinom]{log_lik_multinom()}}, \code{\link[=log_lik_neg_binom]{log_lik_neg_binom()}}, \code{\link[=log_lik_norm]{log_lik_norm()}}, \code{\link[=log_lik_pois]{log_lik_pois()}}, diff --git a/man/log_lik_unif.Rd b/man/log_lik_unif.Rd index e9eef960..90b4fd54 100644 --- a/man/log_lik_unif.Rd +++ b/man/log_lik_unif.Rd @@ -33,6 +33,7 @@ Other log_lik_dist: \code{\link[=log_lik_gamma_pois]{log_lik_gamma_pois()}}, \code{\link[=log_lik_gamma_pois_zi]{log_lik_gamma_pois_zi()}}, \code{\link[=log_lik_lnorm]{log_lik_lnorm()}}, +\code{\link[=log_lik_multinom]{log_lik_multinom()}}, \code{\link[=log_lik_neg_binom]{log_lik_neg_binom()}}, \code{\link[=log_lik_norm]{log_lik_norm()}}, \code{\link[=log_lik_pois]{log_lik_pois()}}, diff --git a/man/params.Rd b/man/params.Rd index 47b31cdf..360f244e 100644 --- a/man/params.Rd +++ b/man/params.Rd @@ -19,6 +19,11 @@ level.} \item{directional}{A flag specifying whether probabilities less than 0.5 should be returned as negative values.} +\item{group}{A vector identifying which rows belong to the same +multinomial trial, i.e., whose \code{x} values must sum to \code{size} and whose +\code{prob} values must sum to 1. Required because the categories within a +single multinomial trial are not independent.} + \item{lambda}{A non-negative numeric vector of means.} \item{level}{A number > 0 and <= 1 specifying the probability coverage of the diff --git a/man/ran_bern.Rd b/man/ran_bern.Rd index bde299dd..20e1613b 100644 --- a/man/ran_bern.Rd +++ b/man/ran_bern.Rd @@ -30,6 +30,7 @@ Other ran_dist: \code{\link[=ran_gamma_pois]{ran_gamma_pois()}}, \code{\link[=ran_gamma_pois_zi]{ran_gamma_pois_zi()}}, \code{\link[=ran_lnorm]{ran_lnorm()}}, +\code{\link[=ran_multinom]{ran_multinom()}}, \code{\link[=ran_neg_binom]{ran_neg_binom()}}, \code{\link[=ran_norm]{ran_norm()}}, \code{\link[=ran_pois]{ran_pois()}}, diff --git a/man/ran_beta_binom.Rd b/man/ran_beta_binom.Rd index daa14527..da22edbb 100644 --- a/man/ran_beta_binom.Rd +++ b/man/ran_beta_binom.Rd @@ -43,6 +43,7 @@ Other ran_dist: \code{\link[=ran_gamma_pois]{ran_gamma_pois()}}, \code{\link[=ran_gamma_pois_zi]{ran_gamma_pois_zi()}}, \code{\link[=ran_lnorm]{ran_lnorm()}}, +\code{\link[=ran_multinom]{ran_multinom()}}, \code{\link[=ran_neg_binom]{ran_neg_binom()}}, \code{\link[=ran_norm]{ran_norm()}}, \code{\link[=ran_pois]{ran_pois()}}, diff --git a/man/ran_binom.Rd b/man/ran_binom.Rd index 3ec3d020..a352e434 100644 --- a/man/ran_binom.Rd +++ b/man/ran_binom.Rd @@ -32,6 +32,7 @@ Other ran_dist: \code{\link[=ran_gamma_pois]{ran_gamma_pois()}}, \code{\link[=ran_gamma_pois_zi]{ran_gamma_pois_zi()}}, \code{\link[=ran_lnorm]{ran_lnorm()}}, +\code{\link[=ran_multinom]{ran_multinom()}}, \code{\link[=ran_neg_binom]{ran_neg_binom()}}, \code{\link[=ran_norm]{ran_norm()}}, \code{\link[=ran_pois]{ran_pois()}}, diff --git a/man/ran_gamma.Rd b/man/ran_gamma.Rd index c348912a..9dc857c7 100644 --- a/man/ran_gamma.Rd +++ b/man/ran_gamma.Rd @@ -31,6 +31,7 @@ Other ran_dist: \code{\link[=ran_gamma_pois]{ran_gamma_pois()}}, \code{\link[=ran_gamma_pois_zi]{ran_gamma_pois_zi()}}, \code{\link[=ran_lnorm]{ran_lnorm()}}, +\code{\link[=ran_multinom]{ran_multinom()}}, \code{\link[=ran_neg_binom]{ran_neg_binom()}}, \code{\link[=ran_norm]{ran_norm()}}, \code{\link[=ran_pois]{ran_pois()}}, diff --git a/man/ran_gamma_pois.Rd b/man/ran_gamma_pois.Rd index 83925a90..342f4c15 100644 --- a/man/ran_gamma_pois.Rd +++ b/man/ran_gamma_pois.Rd @@ -32,6 +32,7 @@ Other ran_dist: \code{\link[=ran_gamma]{ran_gamma()}}, \code{\link[=ran_gamma_pois_zi]{ran_gamma_pois_zi()}}, \code{\link[=ran_lnorm]{ran_lnorm()}}, +\code{\link[=ran_multinom]{ran_multinom()}}, \code{\link[=ran_neg_binom]{ran_neg_binom()}}, \code{\link[=ran_norm]{ran_norm()}}, \code{\link[=ran_pois]{ran_pois()}}, diff --git a/man/ran_gamma_pois_zi.Rd b/man/ran_gamma_pois_zi.Rd index 44240a73..ba438271 100644 --- a/man/ran_gamma_pois_zi.Rd +++ b/man/ran_gamma_pois_zi.Rd @@ -35,6 +35,7 @@ Other ran_dist: \code{\link[=ran_gamma]{ran_gamma()}}, \code{\link[=ran_gamma_pois]{ran_gamma_pois()}}, \code{\link[=ran_lnorm]{ran_lnorm()}}, +\code{\link[=ran_multinom]{ran_multinom()}}, \code{\link[=ran_neg_binom]{ran_neg_binom()}}, \code{\link[=ran_norm]{ran_norm()}}, \code{\link[=ran_pois]{ran_pois()}}, diff --git a/man/ran_lnorm.Rd b/man/ran_lnorm.Rd index e852d444..58b7bb2e 100644 --- a/man/ran_lnorm.Rd +++ b/man/ran_lnorm.Rd @@ -32,6 +32,7 @@ Other ran_dist: \code{\link[=ran_gamma]{ran_gamma()}}, \code{\link[=ran_gamma_pois]{ran_gamma_pois()}}, \code{\link[=ran_gamma_pois_zi]{ran_gamma_pois_zi()}}, +\code{\link[=ran_multinom]{ran_multinom()}}, \code{\link[=ran_neg_binom]{ran_neg_binom()}}, \code{\link[=ran_norm]{ran_norm()}}, \code{\link[=ran_pois]{ran_pois()}}, diff --git a/man/ran_multinom.Rd b/man/ran_multinom.Rd new file mode 100644 index 00000000..bc341737 --- /dev/null +++ b/man/ran_multinom.Rd @@ -0,0 +1,57 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/ran.R +\name{ran_multinom} +\alias{ran_multinom} +\title{Multinomial Random Samples} +\usage{ +ran_multinom(size = 1, prob, group) +} +\arguments{ +\item{size}{A non-negative whole numeric vector of the number of trials.} + +\item{prob}{A numeric vector of the probability of the category. Must sum +to 1 across the rows sharing the same \code{group}.} + +\item{group}{A vector identifying which rows belong to the same +multinomial trial, i.e., whose \code{x} values must sum to \code{size} and whose +\code{prob} values must sum to 1. Required because the categories within a +single multinomial trial are not independent.} +} +\value{ +An integer vector of the random samples, one per row of \code{prob}. +} +\description{ +The multinomial distribution models the counts across two or more +mutually exclusive categories arising from a fixed number of trials. Data +(and therefore random samples) are in \emph{long} format: one value per +category per trial, with \code{group} identifying which rows belong to the +same trial. All rows sharing a \code{group} must have the same \code{size}, and +their \code{prob} values must sum to 1. +} +\details{ +Unlike the other \verb{ran_*()} functions, \code{ran_multinom()} has no \code{n} +argument: the number of samples is fully determined by \code{length(prob)} +(equivalently \code{length(group)}), because a trial's categories can't be +generated independently of one another. +} +\examples{ +ran_multinom(size = 10, prob = c(0.2, 0.3, 0.5), group = c(1, 1, 1)) +} +\seealso{ +Other ran_dist: +\code{\link[=ran_bern]{ran_bern()}}, +\code{\link[=ran_beta_binom]{ran_beta_binom()}}, +\code{\link[=ran_binom]{ran_binom()}}, +\code{\link[=ran_gamma]{ran_gamma()}}, +\code{\link[=ran_gamma_pois]{ran_gamma_pois()}}, +\code{\link[=ran_gamma_pois_zi]{ran_gamma_pois_zi()}}, +\code{\link[=ran_lnorm]{ran_lnorm()}}, +\code{\link[=ran_neg_binom]{ran_neg_binom()}}, +\code{\link[=ran_norm]{ran_norm()}}, +\code{\link[=ran_pois]{ran_pois()}}, +\code{\link[=ran_pois_zi]{ran_pois_zi()}}, +\code{\link[=ran_skewlnorm]{ran_skewlnorm()}}, +\code{\link[=ran_skewnorm]{ran_skewnorm()}}, +\code{\link[=ran_student]{ran_student()}} +} +\concept{ran_dist} diff --git a/man/ran_neg_binom.Rd b/man/ran_neg_binom.Rd index b70d571e..5c66caa5 100644 --- a/man/ran_neg_binom.Rd +++ b/man/ran_neg_binom.Rd @@ -33,6 +33,7 @@ Other ran_dist: \code{\link[=ran_gamma_pois]{ran_gamma_pois()}}, \code{\link[=ran_gamma_pois_zi]{ran_gamma_pois_zi()}}, \code{\link[=ran_lnorm]{ran_lnorm()}}, +\code{\link[=ran_multinom]{ran_multinom()}}, \code{\link[=ran_norm]{ran_norm()}}, \code{\link[=ran_pois]{ran_pois()}}, \code{\link[=ran_pois_zi]{ran_pois_zi()}}, diff --git a/man/ran_norm.Rd b/man/ran_norm.Rd index 84564c6b..1b92e336 100644 --- a/man/ran_norm.Rd +++ b/man/ran_norm.Rd @@ -32,6 +32,7 @@ Other ran_dist: \code{\link[=ran_gamma_pois]{ran_gamma_pois()}}, \code{\link[=ran_gamma_pois_zi]{ran_gamma_pois_zi()}}, \code{\link[=ran_lnorm]{ran_lnorm()}}, +\code{\link[=ran_multinom]{ran_multinom()}}, \code{\link[=ran_neg_binom]{ran_neg_binom()}}, \code{\link[=ran_pois]{ran_pois()}}, \code{\link[=ran_pois_zi]{ran_pois_zi()}}, diff --git a/man/ran_pois.Rd b/man/ran_pois.Rd index 3afe4a9c..e22a041d 100644 --- a/man/ran_pois.Rd +++ b/man/ran_pois.Rd @@ -30,6 +30,7 @@ Other ran_dist: \code{\link[=ran_gamma_pois]{ran_gamma_pois()}}, \code{\link[=ran_gamma_pois_zi]{ran_gamma_pois_zi()}}, \code{\link[=ran_lnorm]{ran_lnorm()}}, +\code{\link[=ran_multinom]{ran_multinom()}}, \code{\link[=ran_neg_binom]{ran_neg_binom()}}, \code{\link[=ran_norm]{ran_norm()}}, \code{\link[=ran_pois_zi]{ran_pois_zi()}}, diff --git a/man/ran_pois_zi.Rd b/man/ran_pois_zi.Rd index 9668cb61..509f5fc0 100644 --- a/man/ran_pois_zi.Rd +++ b/man/ran_pois_zi.Rd @@ -33,6 +33,7 @@ Other ran_dist: \code{\link[=ran_gamma_pois]{ran_gamma_pois()}}, \code{\link[=ran_gamma_pois_zi]{ran_gamma_pois_zi()}}, \code{\link[=ran_lnorm]{ran_lnorm()}}, +\code{\link[=ran_multinom]{ran_multinom()}}, \code{\link[=ran_neg_binom]{ran_neg_binom()}}, \code{\link[=ran_norm]{ran_norm()}}, \code{\link[=ran_pois]{ran_pois()}}, diff --git a/man/ran_skewlnorm.Rd b/man/ran_skewlnorm.Rd index 36f925d1..c02b39d4 100644 --- a/man/ran_skewlnorm.Rd +++ b/man/ran_skewlnorm.Rd @@ -39,6 +39,7 @@ Other ran_dist: \code{\link[=ran_gamma_pois]{ran_gamma_pois()}}, \code{\link[=ran_gamma_pois_zi]{ran_gamma_pois_zi()}}, \code{\link[=ran_lnorm]{ran_lnorm()}}, +\code{\link[=ran_multinom]{ran_multinom()}}, \code{\link[=ran_neg_binom]{ran_neg_binom()}}, \code{\link[=ran_norm]{ran_norm()}}, \code{\link[=ran_pois]{ran_pois()}}, diff --git a/man/ran_skewnorm.Rd b/man/ran_skewnorm.Rd index e8152d03..9b1a4375 100644 --- a/man/ran_skewnorm.Rd +++ b/man/ran_skewnorm.Rd @@ -38,6 +38,7 @@ Other ran_dist: \code{\link[=ran_gamma_pois]{ran_gamma_pois()}}, \code{\link[=ran_gamma_pois_zi]{ran_gamma_pois_zi()}}, \code{\link[=ran_lnorm]{ran_lnorm()}}, +\code{\link[=ran_multinom]{ran_multinom()}}, \code{\link[=ran_neg_binom]{ran_neg_binom()}}, \code{\link[=ran_norm]{ran_norm()}}, \code{\link[=ran_pois]{ran_pois()}}, diff --git a/man/ran_student.Rd b/man/ran_student.Rd index 8a883907..77fdde4c 100644 --- a/man/ran_student.Rd +++ b/man/ran_student.Rd @@ -35,6 +35,7 @@ Other ran_dist: \code{\link[=ran_gamma_pois]{ran_gamma_pois()}}, \code{\link[=ran_gamma_pois_zi]{ran_gamma_pois_zi()}}, \code{\link[=ran_lnorm]{ran_lnorm()}}, +\code{\link[=ran_multinom]{ran_multinom()}}, \code{\link[=ran_neg_binom]{ran_neg_binom()}}, \code{\link[=ran_norm]{ran_norm()}}, \code{\link[=ran_pois]{ran_pois()}}, diff --git a/man/res_bern.Rd b/man/res_bern.Rd index c5beb276..5da2d149 100644 --- a/man/res_bern.Rd +++ b/man/res_bern.Rd @@ -34,6 +34,7 @@ Other res_dist: \code{\link[=res_gamma_pois]{res_gamma_pois()}}, \code{\link[=res_gamma_pois_zi]{res_gamma_pois_zi()}}, \code{\link[=res_lnorm]{res_lnorm()}}, +\code{\link[=res_multinom]{res_multinom()}}, \code{\link[=res_neg_binom]{res_neg_binom()}}, \code{\link[=res_norm]{res_norm()}}, \code{\link[=res_pois]{res_pois()}}, diff --git a/man/res_beta_binom.Rd b/man/res_beta_binom.Rd index 479167cd..5452b575 100644 --- a/man/res_beta_binom.Rd +++ b/man/res_beta_binom.Rd @@ -54,6 +54,7 @@ Other res_dist: \code{\link[=res_gamma_pois]{res_gamma_pois()}}, \code{\link[=res_gamma_pois_zi]{res_gamma_pois_zi()}}, \code{\link[=res_lnorm]{res_lnorm()}}, +\code{\link[=res_multinom]{res_multinom()}}, \code{\link[=res_neg_binom]{res_neg_binom()}}, \code{\link[=res_norm]{res_norm()}}, \code{\link[=res_pois]{res_pois()}}, diff --git a/man/res_binom.Rd b/man/res_binom.Rd index b8026c4e..a00ed641 100644 --- a/man/res_binom.Rd +++ b/man/res_binom.Rd @@ -36,6 +36,7 @@ Other res_dist: \code{\link[=res_gamma_pois]{res_gamma_pois()}}, \code{\link[=res_gamma_pois_zi]{res_gamma_pois_zi()}}, \code{\link[=res_lnorm]{res_lnorm()}}, +\code{\link[=res_multinom]{res_multinom()}}, \code{\link[=res_neg_binom]{res_neg_binom()}}, \code{\link[=res_norm]{res_norm()}}, \code{\link[=res_pois]{res_pois()}}, diff --git a/man/res_gamma.Rd b/man/res_gamma.Rd index 58dff315..6902641e 100644 --- a/man/res_gamma.Rd +++ b/man/res_gamma.Rd @@ -35,6 +35,7 @@ Other res_dist: \code{\link[=res_gamma_pois]{res_gamma_pois()}}, \code{\link[=res_gamma_pois_zi]{res_gamma_pois_zi()}}, \code{\link[=res_lnorm]{res_lnorm()}}, +\code{\link[=res_multinom]{res_multinom()}}, \code{\link[=res_neg_binom]{res_neg_binom()}}, \code{\link[=res_norm]{res_norm()}}, \code{\link[=res_pois]{res_pois()}}, diff --git a/man/res_gamma_pois.Rd b/man/res_gamma_pois.Rd index ae824418..c2e7e481 100644 --- a/man/res_gamma_pois.Rd +++ b/man/res_gamma_pois.Rd @@ -36,6 +36,7 @@ Other res_dist: \code{\link[=res_gamma]{res_gamma()}}, \code{\link[=res_gamma_pois_zi]{res_gamma_pois_zi()}}, \code{\link[=res_lnorm]{res_lnorm()}}, +\code{\link[=res_multinom]{res_multinom()}}, \code{\link[=res_neg_binom]{res_neg_binom()}}, \code{\link[=res_norm]{res_norm()}}, \code{\link[=res_pois]{res_pois()}}, diff --git a/man/res_gamma_pois_zi.Rd b/man/res_gamma_pois_zi.Rd index 5c3080bd..3d4a80d2 100644 --- a/man/res_gamma_pois_zi.Rd +++ b/man/res_gamma_pois_zi.Rd @@ -45,6 +45,7 @@ Other res_dist: \code{\link[=res_gamma]{res_gamma()}}, \code{\link[=res_gamma_pois]{res_gamma_pois()}}, \code{\link[=res_lnorm]{res_lnorm()}}, +\code{\link[=res_multinom]{res_multinom()}}, \code{\link[=res_neg_binom]{res_neg_binom()}}, \code{\link[=res_norm]{res_norm()}}, \code{\link[=res_pois]{res_pois()}}, diff --git a/man/res_lnorm.Rd b/man/res_lnorm.Rd index 3fbfd1a5..1a674850 100644 --- a/man/res_lnorm.Rd +++ b/man/res_lnorm.Rd @@ -36,6 +36,7 @@ Other res_dist: \code{\link[=res_gamma]{res_gamma()}}, \code{\link[=res_gamma_pois]{res_gamma_pois()}}, \code{\link[=res_gamma_pois_zi]{res_gamma_pois_zi()}}, +\code{\link[=res_multinom]{res_multinom()}}, \code{\link[=res_neg_binom]{res_neg_binom()}}, \code{\link[=res_norm]{res_norm()}}, \code{\link[=res_pois]{res_pois()}}, diff --git a/man/res_multinom.Rd b/man/res_multinom.Rd new file mode 100644 index 00000000..53653b68 --- /dev/null +++ b/man/res_multinom.Rd @@ -0,0 +1,77 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/res.R +\name{res_multinom} +\alias{res_multinom} +\title{Multinomial Residuals} +\usage{ +res_multinom(x, size = 1, prob, group, type = "dev", simulate = FALSE) +} +\arguments{ +\item{x}{A non-negative whole numeric vector of the category counts.} + +\item{size}{A non-negative whole numeric vector of the number of trials.} + +\item{prob}{A numeric vector of the probability of the category. Must sum +to 1 across the rows sharing the same \code{group}.} + +\item{group}{A vector identifying which rows belong to the same +multinomial trial, i.e., whose \code{x} values must sum to \code{size} and whose +\code{prob} values must sum to 1. Required because the categories within a +single multinomial trial are not independent.} + +\item{type}{A string of the residual type. 'raw' for raw residuals 'dev' for +deviance residuals and 'data' for the data.} + +\item{simulate}{A flag specifying whether to simulate residuals.} +} +\value{ +An numeric vector of the corresponding residuals. +} +\description{ +The multinomial distribution models the counts across two or more +mutually exclusive categories arising from a fixed number of trials. Data +are in \emph{long} format: one row per category per trial, with \code{group} +identifying which rows belong to the same trial (see +\code{\link[=log_lik_multinom]{log_lik_multinom()}} for details). \code{res_multinom()} returns one residual +per row (per category per trial), not one per trial -- since a trial's +categories are not independent, there's no single meaningful residual (or +sign) for the trial as a whole. The classic per-trial deviance statistic +can always be recovered by summing the squared \code{type = "dev"} residuals +within a \code{group}. +} +\details{ +Because a category count is marginally binomial (\code{x ~ Binomial(size, prob)}), the \code{"standardized"} residual uses the same formula as +\code{\link[=res_binom]{res_binom()}}; this is already the covariance-adjusted (unit-variance) +residual, not the naive independent-Poisson version. + +\code{group} is only used when \code{simulate = TRUE}, to draw a joint, +correlation-preserving replicate for each trial (via \code{\link[=ran_multinom]{ran_multinom()}}) +rather than simulating each category independently, which would +understate the true (negative) covariance among a trial's categories. For +the simulated residuals from this function to be valid inputs to +\code{embr::posterior_predictive_check()}, the model must be fit with +\code{new_expr_vec = TRUE}, so that \code{res_multinom()} is evaluated once on the +full data vector (with visibility of every row in each \code{group}) rather +than row by row. +} +\examples{ +res_multinom(c(1, 3, 6), size = 10, prob = c(0.2, 0.3, 0.5), group = c(1, 1, 1)) +} +\seealso{ +Other res_dist: +\code{\link[=res_bern]{res_bern()}}, +\code{\link[=res_beta_binom]{res_beta_binom()}}, +\code{\link[=res_binom]{res_binom()}}, +\code{\link[=res_gamma]{res_gamma()}}, +\code{\link[=res_gamma_pois]{res_gamma_pois()}}, +\code{\link[=res_gamma_pois_zi]{res_gamma_pois_zi()}}, +\code{\link[=res_lnorm]{res_lnorm()}}, +\code{\link[=res_neg_binom]{res_neg_binom()}}, +\code{\link[=res_norm]{res_norm()}}, +\code{\link[=res_pois]{res_pois()}}, +\code{\link[=res_pois_zi]{res_pois_zi()}}, +\code{\link[=res_skewlnorm]{res_skewlnorm()}}, +\code{\link[=res_skewnorm]{res_skewnorm()}}, +\code{\link[=res_student]{res_student()}} +} +\concept{res_dist} diff --git a/man/res_neg_binom.Rd b/man/res_neg_binom.Rd index 50cf1915..f92338ba 100644 --- a/man/res_neg_binom.Rd +++ b/man/res_neg_binom.Rd @@ -37,6 +37,7 @@ Other res_dist: \code{\link[=res_gamma_pois]{res_gamma_pois()}}, \code{\link[=res_gamma_pois_zi]{res_gamma_pois_zi()}}, \code{\link[=res_lnorm]{res_lnorm()}}, +\code{\link[=res_multinom]{res_multinom()}}, \code{\link[=res_norm]{res_norm()}}, \code{\link[=res_pois]{res_pois()}}, \code{\link[=res_pois_zi]{res_pois_zi()}}, diff --git a/man/res_norm.Rd b/man/res_norm.Rd index ebabf927..3c6eff54 100644 --- a/man/res_norm.Rd +++ b/man/res_norm.Rd @@ -36,6 +36,7 @@ Other res_dist: \code{\link[=res_gamma_pois]{res_gamma_pois()}}, \code{\link[=res_gamma_pois_zi]{res_gamma_pois_zi()}}, \code{\link[=res_lnorm]{res_lnorm()}}, +\code{\link[=res_multinom]{res_multinom()}}, \code{\link[=res_neg_binom]{res_neg_binom()}}, \code{\link[=res_pois]{res_pois()}}, \code{\link[=res_pois_zi]{res_pois_zi()}}, diff --git a/man/res_pois.Rd b/man/res_pois.Rd index 6da9d015..193bbea0 100644 --- a/man/res_pois.Rd +++ b/man/res_pois.Rd @@ -34,6 +34,7 @@ Other res_dist: \code{\link[=res_gamma_pois]{res_gamma_pois()}}, \code{\link[=res_gamma_pois_zi]{res_gamma_pois_zi()}}, \code{\link[=res_lnorm]{res_lnorm()}}, +\code{\link[=res_multinom]{res_multinom()}}, \code{\link[=res_neg_binom]{res_neg_binom()}}, \code{\link[=res_norm]{res_norm()}}, \code{\link[=res_pois_zi]{res_pois_zi()}}, diff --git a/man/res_pois_zi.Rd b/man/res_pois_zi.Rd index 8b2aaea2..c8706f20 100644 --- a/man/res_pois_zi.Rd +++ b/man/res_pois_zi.Rd @@ -36,6 +36,7 @@ Other res_dist: \code{\link[=res_gamma_pois]{res_gamma_pois()}}, \code{\link[=res_gamma_pois_zi]{res_gamma_pois_zi()}}, \code{\link[=res_lnorm]{res_lnorm()}}, +\code{\link[=res_multinom]{res_multinom()}}, \code{\link[=res_neg_binom]{res_neg_binom()}}, \code{\link[=res_norm]{res_norm()}}, \code{\link[=res_pois]{res_pois()}}, diff --git a/man/res_skewlnorm.Rd b/man/res_skewlnorm.Rd index 98364d56..535b3bee 100644 --- a/man/res_skewlnorm.Rd +++ b/man/res_skewlnorm.Rd @@ -48,6 +48,7 @@ Other res_dist: \code{\link[=res_gamma_pois]{res_gamma_pois()}}, \code{\link[=res_gamma_pois_zi]{res_gamma_pois_zi()}}, \code{\link[=res_lnorm]{res_lnorm()}}, +\code{\link[=res_multinom]{res_multinom()}}, \code{\link[=res_neg_binom]{res_neg_binom()}}, \code{\link[=res_norm]{res_norm()}}, \code{\link[=res_pois]{res_pois()}}, diff --git a/man/res_skewnorm.Rd b/man/res_skewnorm.Rd index f9ec34cf..7af499f1 100644 --- a/man/res_skewnorm.Rd +++ b/man/res_skewnorm.Rd @@ -40,6 +40,7 @@ Other res_dist: \code{\link[=res_gamma_pois]{res_gamma_pois()}}, \code{\link[=res_gamma_pois_zi]{res_gamma_pois_zi()}}, \code{\link[=res_lnorm]{res_lnorm()}}, +\code{\link[=res_multinom]{res_multinom()}}, \code{\link[=res_neg_binom]{res_neg_binom()}}, \code{\link[=res_norm]{res_norm()}}, \code{\link[=res_pois]{res_pois()}}, diff --git a/man/res_student.Rd b/man/res_student.Rd index 9f9bf7b8..13d5c8cd 100644 --- a/man/res_student.Rd +++ b/man/res_student.Rd @@ -39,6 +39,7 @@ Other res_dist: \code{\link[=res_gamma_pois]{res_gamma_pois()}}, \code{\link[=res_gamma_pois_zi]{res_gamma_pois_zi()}}, \code{\link[=res_lnorm]{res_lnorm()}}, +\code{\link[=res_multinom]{res_multinom()}}, \code{\link[=res_neg_binom]{res_neg_binom()}}, \code{\link[=res_norm]{res_norm()}}, \code{\link[=res_pois]{res_pois()}}, diff --git a/tests/testthat/test-dev.R b/tests/testthat/test-dev.R index 1afb5c5e..8351d1a8 100644 --- a/tests/testthat/test-dev.R +++ b/tests/testthat/test-dev.R @@ -473,6 +473,79 @@ test_that("dev_lnorm deviance", { expect_equal(deviance, deviance(mod)) }) +test_that("dev_multinom", { + expect_identical( + dev_multinom(numeric(0), numeric(0), numeric(0)), + numeric(0) + ) + expect_identical( + dev_multinom(c(2, 3, 5), size = 10, prob = c(0.2, 0.3, 0.5)), + c(0, 0, 0) + ) + expect_identical( + dev_multinom(1, 10, 0.2), + dev_pois(1, 10 * 0.2) + ) + # row-level deviance identity: summing the per-cell deviances within a + # trial recovers the classic multinomial saturated-model deviance + x <- c(1, 3, 6) + size <- 10 + prob <- c(0.2, 0.3, 0.5) + expect_equal( + sum(dev_multinom(x, size, prob)), + 2 * sum(x * log(x / (size * prob))) + ) +}) + +test_that("dev_multinom matches glmnet deviance", { + skip_if_not_installed("glmnet") + withr::with_seed(7, { + n <- 150 + K <- 4 + x <- cbind(rnorm(n), rnorm(n), rnorm(n)) + eta <- cbind( + 0, + 0.9 * x[, 1] - 0.4 * x[, 2], + -0.3 * x[, 1] + 0.6 * x[, 3], + 0.5 * x[, 2] - 0.2 * x[, 3] + ) + prob_true <- exp(eta) / rowSums(exp(eta)) + y_wide <- t(sapply(seq_len(n), function(i) { + stats::rmultinom(1, size = 1, prob = prob_true[i, ]) + })) + colnames(y_wide) <- paste0("cat", seq_len(K)) + + for (lambda in c(0, 0.01, 0.05)) { + fit <- glmnet::glmnet( + x, + y_wide, + family = "multinomial", + lambda = lambda, + intercept = TRUE + ) + dev_glmnet <- deviance(fit) + prob_hat <- predict(fit, newx = x, type = "response")[, , 1] + expect_equal( + sum(dev_multinom( + as.vector(t(y_wide)), + size = 1, + prob = as.vector(t(prob_hat)) + )), + as.numeric(dev_glmnet) + ) + } + }) +}) + +test_that("dev_multinom res", { + x <- c(1, 3, 6) + size <- 10 + prob <- c(0.2, 0.3, 0.5) + res <- dev_multinom(x, size, prob, res = TRUE) + expect_equal(sign(res), sign(x - size * prob)) + expect_equal(sum(res^2), sum(dev_multinom(x, size, prob))) +}) + test_that("dev_neg_binom", { expect_identical( dev_neg_binom(integer(0), integer(0), integer(0)), diff --git a/tests/testthat/test-log-lik.R b/tests/testthat/test-log-lik.R index 863441ca..99841b87 100644 --- a/tests/testthat/test-log-lik.R +++ b/tests/testthat/test-log-lik.R @@ -332,6 +332,44 @@ test_that("log_lik_lnorm truncated", { ) }) +test_that("log_lik_multinom", { + expect_identical( + log_lik_multinom(numeric(0), numeric(0), numeric(0), numeric(0)), + numeric(0) + ) + expect_error( + log_lik_multinom(c(1, 3, 6), c(10, 10, 5), c(0.2, 0.3, 0.5), c(1, 1, 1)), + "`size` must be the same for every row belonging to the same `group`" + ) + expect_error( + log_lik_multinom(c(1, 3, 6), 10, c(0.2, 0.3, 0.4), c(1, 1, 1)), + "`prob` must sum to 1 for every `group`" + ) + expect_equal( + sum(log_lik_multinom( + c(1, 3, 6), + size = 10, + prob = c(0.2, 0.3, 0.5), + group = c(1, 1, 1) + )), + dmultinom(c(1, 3, 6), size = 10, prob = c(0.2, 0.3, 0.5), log = TRUE) + ) + # multiple trials in long format + x <- c(1, 3, 6, 2, 2) + size <- c(10, 10, 10, 4, 4) + prob <- c(0.2, 0.3, 0.5, 0.5, 0.5) + group <- c(1, 1, 1, 2, 2) + ll <- log_lik_multinom(x, size, prob, group) + expect_equal( + sum(ll[group == 1]), + dmultinom(x[group == 1], size = 10, prob = prob[group == 1], log = TRUE) + ) + expect_equal( + sum(ll[group == 2]), + dmultinom(x[group == 2], size = 4, prob = prob[group == 2], log = TRUE) + ) +}) + test_that("log_lik_neg_binom", { expect_identical( log_lik_neg_binom(0, 2, 1), diff --git a/tests/testthat/test-ran.R b/tests/testthat/test-ran.R index 4fcc5e00..d9882fe1 100644 --- a/tests/testthat/test-ran.R +++ b/tests/testthat/test-ran.R @@ -107,6 +107,35 @@ test_that("ran_gamma_pois", { }) }) +test_that("ran_multinom", { + expect_identical( + ran_multinom(size = numeric(0), prob = numeric(0), group = numeric(0)), + integer(0) + ) + expect_error( + ran_multinom(size = c(10, 5), prob = c(0.5, 0.5), group = c(1, 1)), + "`size` must be the same for every row belonging to the same `group`" + ) + expect_error( + ran_multinom(size = 10, prob = c(0.5, 0.4), group = c(1, 1)), + "`prob` must sum to 1 for every `group`" + ) + withr::with_seed(101, { + x <- ran_multinom(size = 10, prob = c(0.2, 0.3, 0.5), group = c(1, 1, 1)) + expect_identical(sum(x), 10L) + expect_length(x, 3L) + }) + withr::with_seed(101, { + x <- ran_multinom( + size = c(10, 10, 6, 6), + prob = c(0.2, 0.8, 0.5, 0.5), + group = c(1, 1, 2, 2) + ) + expect_identical(x[1] + x[2], 10L) + expect_identical(x[3] + x[4], 6L) + }) +}) + test_that("ran_neg_binom", { expect_error(ran_neg_binom(NA_integer_)) expect_error(ran_neg_binom(integer(0))) diff --git a/tests/testthat/test-res.R b/tests/testthat/test-res.R index e7506af7..8145acb0 100644 --- a/tests/testthat/test-res.R +++ b/tests/testthat/test-res.R @@ -507,6 +507,45 @@ test_that("res_gamma_pois simulate", { }) }) +test_that("res_multinom", { + x <- c(1, 3, 6) + size <- 10 + prob <- c(0.2, 0.3, 0.5) + group <- c(1, 1, 1) + + expect_equal( + res_multinom(x, size, prob, group, type = "dev"), + dev_multinom(x, size, prob, res = TRUE) + ) + expect_equal( + res_multinom(x, size, prob, group, type = "raw"), + x - size * prob + ) + # a category count is marginally binomial, so the standardized residual + # matches res_binom exactly + expect_equal( + res_multinom(x, size, prob, group, type = "standardized"), + res_binom(x, size, prob, type = "standardized") + ) + expect_equal( + res_multinom(x, size, prob, group, type = "data"), + x + ) + expect_error(res_multinom(x, size, prob, group, type = "unknown")) + + # sum of squared deviance residuals recovers the row-level deviance + expect_equal( + sum(res_multinom(x, size, prob, group, type = "dev")^2), + 2 * sum(x * log(x / (size * prob))) + ) + + withr::with_seed(101, { + sim <- res_multinom(x, size, prob, group, type = "data", simulate = TRUE) + expect_identical(sum(sim), as.integer(size)) + expect_length(sim, 3L) + }) +}) + test_that("res_neg_binom", { expect_identical( res_neg_binom(integer(0), integer(0), integer(0)), From 2339606c59a13adf3e7b8e6088b43df178b283db Mon Sep 17 00:00:00 2001 From: Nicole Hill Date: Thu, 6 Aug 2026 15:26:07 -0700 Subject: [PATCH 2/9] Fix log_lik_multinom() ave() bug and guard against singleton groups ave()'s non-standard FUN lookup fails when log_lik_multinom() is evaluated inside a Bayesian model-fitting framework's derived-quantity expressions; replaced with an equivalent table()-based computation. Also reject `group` values mapping to a single row in chk_multinom_group() (shared by ran_multinom()/log_lik_multinom(), and so reached by res_multinom(simulate = TRUE)): a singleton group can't represent a real multinomial trial, and silently produces a degenerate result (rmultinom() just returns `size`) rather than erroring. This is the shape you get if the function ends up evaluated once per row instead of once for the whole data vector. --- R/internal.R | 13 ++++++++++++- R/log-lik.R | 3 ++- R/res.R | 10 ++++------ man/res_multinom.Rd | 10 ++++------ tests/testthat/test-log-lik.R | 4 ++++ tests/testthat/test-ran.R | 12 ++++++++++++ 6 files changed, 38 insertions(+), 14 deletions(-) diff --git a/R/internal.R b/R/internal.R index 9257b6a8..008416e1 100644 --- a/R/internal.R +++ b/R/internal.R @@ -5,9 +5,20 @@ dev_res <- function(x, mu, dev) { # Validates that every row sharing a `group` (multinomial trial) has the # same `size` and that its `prob` values sum to 1, since both the # Poisson-equivalent deviance/log-likelihood identities and simulation via -# rmultinom() require this to hold. +# rmultinom() require this to hold. Also rejects singleton groups: a group +# with only one row can't represent a real multinomial trial (a trial needs +# at least 2 categories) -- rmultinom() would just deterministically return +# `size` for it. This typically happens if the function ends up being +# evaluated separately for each row instead of once for the whole data +# vector, so each call only ever sees a single row's own group value. chk_multinom_group <- function(size, prob, group) { for (idx in split(seq_along(group), group)) { + if (length(idx) < 2L) { + stop( + "Each `group` must contain at least 2 rows (a multinomial trial needs at least 2 categories); found a group with only 1 row. This usually means `group` lost the rows that should share a trial before reaching here.", + call. = FALSE + ) + } if (length(unique(size[idx])) > 1L) { stop( "`size` must be the same for every row belonging to the same `group` (multinomial trial).", diff --git a/R/log-lik.R b/R/log-lik.R index eb83b436..6b020810 100644 --- a/R/log-lik.R +++ b/R/log-lik.R @@ -468,7 +468,8 @@ log_lik_multinom <- function(x, size = 1, prob, group) { chk_multinom_group(size, prob, group) mu <- size * prob log_lik <- log_lik_pois(x, mu) - k <- ave(seq_along(group), group, FUN = length) + group_size <- table(group) + k <- as.numeric(group_size[as.character(group)]) const <- log_lik_pois(size, size) log_lik - const / k } diff --git a/R/res.R b/R/res.R index 2b701063..16f71c66 100644 --- a/R/res.R +++ b/R/res.R @@ -272,12 +272,10 @@ res_lnorm <- function( #' `group` is only used when `simulate = TRUE`, to draw a joint, #' correlation-preserving replicate for each trial (via [ran_multinom()]) #' rather than simulating each category independently, which would -#' understate the true (negative) covariance among a trial's categories. For -#' the simulated residuals from this function to be valid inputs to -#' `embr::posterior_predictive_check()`, the model must be fit with -#' `new_expr_vec = TRUE`, so that `res_multinom()` is evaluated once on the -#' full data vector (with visibility of every row in each `group`) rather -#' than row by row. +#' understate the true (negative) covariance among a trial's categories. +#' This requires `res_multinom()` to be evaluated on the full data vector at +#' once, so that every row sharing a `group` is visible in the same call, +#' rather than evaluated separately for each row. #' #' @inheritParams params #' @param x A non-negative whole numeric vector of the category counts. diff --git a/man/res_multinom.Rd b/man/res_multinom.Rd index 53653b68..0d807e30 100644 --- a/man/res_multinom.Rd +++ b/man/res_multinom.Rd @@ -47,12 +47,10 @@ residual, not the naive independent-Poisson version. \code{group} is only used when \code{simulate = TRUE}, to draw a joint, correlation-preserving replicate for each trial (via \code{\link[=ran_multinom]{ran_multinom()}}) rather than simulating each category independently, which would -understate the true (negative) covariance among a trial's categories. For -the simulated residuals from this function to be valid inputs to -\code{embr::posterior_predictive_check()}, the model must be fit with -\code{new_expr_vec = TRUE}, so that \code{res_multinom()} is evaluated once on the -full data vector (with visibility of every row in each \code{group}) rather -than row by row. +understate the true (negative) covariance among a trial's categories. +This requires \code{res_multinom()} to be evaluated on the full data vector at +once, so that every row sharing a \code{group} is visible in the same call, +rather than evaluated separately for each row. } \examples{ res_multinom(c(1, 3, 6), size = 10, prob = c(0.2, 0.3, 0.5), group = c(1, 1, 1)) diff --git a/tests/testthat/test-log-lik.R b/tests/testthat/test-log-lik.R index 99841b87..01549db9 100644 --- a/tests/testthat/test-log-lik.R +++ b/tests/testthat/test-log-lik.R @@ -345,6 +345,10 @@ test_that("log_lik_multinom", { log_lik_multinom(c(1, 3, 6), 10, c(0.2, 0.3, 0.4), c(1, 1, 1)), "`prob` must sum to 1 for every `group`" ) + expect_error( + log_lik_multinom(4, 10, 1, 1), + "must contain at least 2 rows" + ) expect_equal( sum(log_lik_multinom( c(1, 3, 6), diff --git a/tests/testthat/test-ran.R b/tests/testthat/test-ran.R index d9882fe1..be711fcf 100644 --- a/tests/testthat/test-ran.R +++ b/tests/testthat/test-ran.R @@ -120,6 +120,18 @@ test_that("ran_multinom", { ran_multinom(size = 10, prob = c(0.5, 0.4), group = c(1, 1)), "`prob` must sum to 1 for every `group`" ) + expect_error( + ran_multinom(size = 10, prob = 1, group = 1), + "must contain at least 2 rows" + ) + expect_error( + ran_multinom( + size = c(10, 10, 10, 10), + prob = c(0.2, 0.3, 0.5, 1), + group = c(1, 1, 1, 2) + ), + "must contain at least 2 rows" + ) withr::with_seed(101, { x <- ran_multinom(size = 10, prob = c(0.2, 0.3, 0.5), group = c(1, 1, 1)) expect_identical(sum(x), 10L) From 209243ffaa78a96dfa287a923baf555ed04dcae7 Mon Sep 17 00:00:00 2001 From: Nicole Hill Date: Thu, 6 Aug 2026 15:40:59 -0700 Subject: [PATCH 3/9] Document group's min-rows requirement; add res_multinom() simulate test Note in the shared group param docs that each group must contain at least 2 rows, matching the check already enforced in chk_multinom_group(). Also drop an unneeded implementation aside from res_multinom()'s docs, and add a simulate-many-and-check-mean/sd test for res_multinom(), matching the pattern used for every other distribution's res_*() function. --- R/params.R | 5 +++-- R/res.R | 5 ----- man/log_lik_multinom.Rd | 5 +++-- man/params.Rd | 5 +++-- man/ran_multinom.Rd | 5 +++-- man/res_multinom.Rd | 9 +++------ tests/testthat/test-res.R | 32 ++++++++++++++++++++++++++++++++ 7 files changed, 47 insertions(+), 19 deletions(-) diff --git a/R/params.R b/R/params.R index 513fb8a5..502cb488 100644 --- a/R/params.R +++ b/R/params.R @@ -18,8 +18,9 @@ #' should be returned as negative values. #' @param group A vector identifying which rows belong to the same #' multinomial trial, i.e., whose `x` values must sum to `size` and whose -#' `prob` values must sum to 1. Required because the categories within a -#' single multinomial trial are not independent. +#' `prob` values must sum to 1. Each group must contain at least 2 rows (a +#' trial needs at least 2 categories). Required because the categories +#' within a single multinomial trial are not independent. #' @param lambda A non-negative numeric vector of means. #' @param level A number > 0 and <= 1 specifying the probability coverage of the #' interval. diff --git a/R/res.R b/R/res.R index 16f71c66..02a46df4 100644 --- a/R/res.R +++ b/R/res.R @@ -264,11 +264,6 @@ res_lnorm <- function( #' can always be recovered by summing the squared `type = "dev"` residuals #' within a `group`. #' -#' Because a category count is marginally binomial (`x ~ Binomial(size, -#' prob)`), the `"standardized"` residual uses the same formula as -#' [res_binom()]; this is already the covariance-adjusted (unit-variance) -#' residual, not the naive independent-Poisson version. -#' #' `group` is only used when `simulate = TRUE`, to draw a joint, #' correlation-preserving replicate for each trial (via [ran_multinom()]) #' rather than simulating each category independently, which would diff --git a/man/log_lik_multinom.Rd b/man/log_lik_multinom.Rd index 65fc845d..7367d3ef 100644 --- a/man/log_lik_multinom.Rd +++ b/man/log_lik_multinom.Rd @@ -16,8 +16,9 @@ to 1 across the rows sharing the same \code{group}.} \item{group}{A vector identifying which rows belong to the same multinomial trial, i.e., whose \code{x} values must sum to \code{size} and whose -\code{prob} values must sum to 1. Required because the categories within a -single multinomial trial are not independent.} +\code{prob} values must sum to 1. Each group must contain at least 2 rows (a +trial needs at least 2 categories). Required because the categories +within a single multinomial trial are not independent.} } \value{ An numeric vector of the corresponding log-likelihoods, one value diff --git a/man/params.Rd b/man/params.Rd index 360f244e..b2f7d92f 100644 --- a/man/params.Rd +++ b/man/params.Rd @@ -21,8 +21,9 @@ should be returned as negative values.} \item{group}{A vector identifying which rows belong to the same multinomial trial, i.e., whose \code{x} values must sum to \code{size} and whose -\code{prob} values must sum to 1. Required because the categories within a -single multinomial trial are not independent.} +\code{prob} values must sum to 1. Each group must contain at least 2 rows (a +trial needs at least 2 categories). Required because the categories +within a single multinomial trial are not independent.} \item{lambda}{A non-negative numeric vector of means.} diff --git a/man/ran_multinom.Rd b/man/ran_multinom.Rd index bc341737..dd0842e7 100644 --- a/man/ran_multinom.Rd +++ b/man/ran_multinom.Rd @@ -14,8 +14,9 @@ to 1 across the rows sharing the same \code{group}.} \item{group}{A vector identifying which rows belong to the same multinomial trial, i.e., whose \code{x} values must sum to \code{size} and whose -\code{prob} values must sum to 1. Required because the categories within a -single multinomial trial are not independent.} +\code{prob} values must sum to 1. Each group must contain at least 2 rows (a +trial needs at least 2 categories). Required because the categories +within a single multinomial trial are not independent.} } \value{ An integer vector of the random samples, one per row of \code{prob}. diff --git a/man/res_multinom.Rd b/man/res_multinom.Rd index 0d807e30..15672878 100644 --- a/man/res_multinom.Rd +++ b/man/res_multinom.Rd @@ -16,8 +16,9 @@ to 1 across the rows sharing the same \code{group}.} \item{group}{A vector identifying which rows belong to the same multinomial trial, i.e., whose \code{x} values must sum to \code{size} and whose -\code{prob} values must sum to 1. Required because the categories within a -single multinomial trial are not independent.} +\code{prob} values must sum to 1. Each group must contain at least 2 rows (a +trial needs at least 2 categories). Required because the categories +within a single multinomial trial are not independent.} \item{type}{A string of the residual type. 'raw' for raw residuals 'dev' for deviance residuals and 'data' for the data.} @@ -40,10 +41,6 @@ can always be recovered by summing the squared \code{type = "dev"} residuals within a \code{group}. } \details{ -Because a category count is marginally binomial (\code{x ~ Binomial(size, prob)}), the \code{"standardized"} residual uses the same formula as -\code{\link[=res_binom]{res_binom()}}; this is already the covariance-adjusted (unit-variance) -residual, not the naive independent-Poisson version. - \code{group} is only used when \code{simulate = TRUE}, to draw a joint, correlation-preserving replicate for each trial (via \code{\link[=ran_multinom]{ran_multinom()}}) rather than simulating each category independently, which would diff --git a/tests/testthat/test-res.R b/tests/testthat/test-res.R index 8145acb0..8a376264 100644 --- a/tests/testthat/test-res.R +++ b/tests/testthat/test-res.R @@ -546,6 +546,38 @@ test_that("res_multinom", { }) }) +test_that("res_multinom simulate", { + n_group <- 10000 + size <- rep(10, n_group * 3) + prob <- rep(c(0.2, 0.3, 0.5), n_group) + group <- rep(seq_len(n_group), each = 3) + + withr::with_seed(101, { + res <- res_multinom( + rep(0, n_group * 3), + size = size, + prob = prob, + group = group, + simulate = TRUE, + type = "dev" + ) + expect_equal(mean(res), -0.0794000903313601) + expect_equal(sd(res), 0.859943455555756) + }) + withr::with_seed(101, { + res <- res_multinom( + rep(0, n_group * 3), + size = size, + prob = prob, + group = group, + simulate = TRUE, + type = "standardized" + ) + expect_equal(mean(res), 0.000255549022088425) + expect_equal(sd(res), 1.00164264126585) + }) +}) + test_that("res_neg_binom", { expect_identical( res_neg_binom(integer(0), integer(0), integer(0)), From 9f33a210fc5f18997a2abf2c145c88faf38a5051 Mon Sep 17 00:00:00 2001 From: Nicole Hill Date: Thu, 6 Aug 2026 15:52:35 -0700 Subject: [PATCH 4/9] Test that res_multinom(simulate = TRUE) propagates the singleton-group error --- tests/testthat/test-res.R | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/tests/testthat/test-res.R b/tests/testthat/test-res.R index 8a376264..715f66f5 100644 --- a/tests/testthat/test-res.R +++ b/tests/testthat/test-res.R @@ -532,6 +532,10 @@ test_that("res_multinom", { x ) expect_error(res_multinom(x, size, prob, group, type = "unknown")) + expect_error( + res_multinom(4, 10, 1, 1, simulate = TRUE), + "must contain at least 2 rows" + ) # sum of squared deviance residuals recovers the row-level deviance expect_equal( From 7917cc3ff2b23760e6032ded04a4c54973d71ec1 Mon Sep 17 00:00:00 2001 From: Nicole Hill Date: Fri, 7 Aug 2026 08:43:46 -0700 Subject: [PATCH 5/9] Guard against uneven group sizes and NA group in multinomial functions Reject groups whose row count differs from the mode across the call (a trial silently missing a row for ordinary multinomial logistic regression), and move the NA-group check into each calling function so it runs before group is used to split the data. --- R/internal.R | 65 +++++++++++++++++++++++++++++++++-- R/log-lik.R | 5 ++- R/params.R | 13 +++++-- R/ran.R | 10 +++++- _pkgdown.yml | 4 +++ man/log_lik_multinom.Rd | 13 +++++-- man/params.Rd | 13 +++++-- man/ran_multinom.Rd | 18 +++++++--- man/res_multinom.Rd | 13 +++++-- tests/testthat/test-log-lik.R | 53 +++++++++++++++++++++++++--- tests/testthat/test-ran.R | 43 +++++++++++++++++++++++ tests/testthat/test-res.R | 43 +++++++++++++++++++++++ 12 files changed, 268 insertions(+), 25 deletions(-) diff --git a/R/internal.R b/R/internal.R index 008416e1..9d3bac5c 100644 --- a/R/internal.R +++ b/R/internal.R @@ -11,15 +11,41 @@ dev_res <- function(x, mu, dev) { # `size` for it. This typically happens if the function ends up being # evaluated separately for each row instead of once for the whole data # vector, so each call only ever sees a single row's own group value. +# +# Callers must reject NA `group` values (chk_not_any_na(group)) before +# calling this: unlike `size`/`prob` (attributes of an already-identified +# trial, where a missing value leaves that one trial undefined -- see +# multinom_row_na()), `group` is what identifies the trial in the first +# place. There's no sensible way to know which trial (if any) an unlabelled +# row belongs to, so it can't be resolved by attributing NA to "the trial it +# belongs to" the way size/prob can. That check lives in each calling +# function rather than here so it runs before `group` is used to split the +# data below. +# +# The size/prob checks compare only the non-NA values within a group, so a +# merely-missing value doesn't trigger a spurious error -- it's handled by +# multinom_row_na() instead, which marks the whole group's output NA rather +# than erroring (see its own comment for why). +# +# Also compares every group's row count (its number of categories) against +# the most common row count across all groups in the call, erroring on any +# that differ. This is specifically for ordinary multinomial logistic +# regression, where every trial has the same fixed set of possible +# categories: a group with fewer rows than the rest of the dataset is +# almost always a trial silently missing a row, which would otherwise be +# validated/simulated/scored as if it were a genuinely smaller, complete +# trial, giving a wrong (not NA, not an error) answer. chk_multinom_group <- function(size, prob, group) { - for (idx in split(seq_along(group), group)) { + groups <- split(seq_along(group), group) + for (idx in groups) { if (length(idx) < 2L) { stop( "Each `group` must contain at least 2 rows (a multinomial trial needs at least 2 categories); found a group with only 1 row. This usually means `group` lost the rows that should share a trial before reaching here.", call. = FALSE ) } - if (length(unique(size[idx])) > 1L) { + known_size <- size[idx][!is.na(size[idx])] + if (length(unique(known_size)) > 1L) { stop( "`size` must be the same for every row belonging to the same `group` (multinomial trial).", call. = FALSE @@ -33,4 +59,39 @@ chk_multinom_group <- function(size, prob, group) { ) } } + if (length(groups) > 1L) { + group_sizes <- lengths(groups) + size_counts <- table(group_sizes) + mode_size <- as.integer(names(size_counts)[which.max(size_counts)]) + bad <- group_sizes != mode_size + if (any(bad)) { + stop( + sprintf( + "Every `group` should have the same number of rows (%d, the most common number of categories in this data); found a group (\"%s\") with %d row(s) instead. This usually means `group` lost a row that should have been part of that trial.", + mode_size, + names(groups)[bad][1], + group_sizes[bad][1] + ), + call. = FALSE + ) + } + } +} + +# Returns, for each row, whether its multinomial trial's draw/score can't be +# computed because `size`/`prob` is NA for ANY row sharing its group. A +# multinomial trial is one joint computation across all of its rows (the +# categories aren't independent), so missing information anywhere in that +# trial leaves the whole trial's result undefined -- not just the row that +# happens to be missing a value. `group` itself can't be NA (rejected by +# chk_multinom_group()), so every row here belongs to a real, known trial. +multinom_row_na <- function(size, prob, group) { + bad <- is.na(size) | is.na(prob) + result <- rep(FALSE, length(group)) + for (idx in split(seq_along(group), group)) { + if (any(bad[idx])) { + result[idx] <- TRUE + } + } + result } diff --git a/R/log-lik.R b/R/log-lik.R index 6b020810..f5c646d5 100644 --- a/R/log-lik.R +++ b/R/log-lik.R @@ -465,13 +465,16 @@ log_lik_multinom <- function(x, size = 1, prob, group) { size <- rep_len(size, n) prob <- rep_len(prob, n) group <- rep_len(group, n) + chk_not_any_na(group) chk_multinom_group(size, prob, group) mu <- size * prob log_lik <- log_lik_pois(x, mu) group_size <- table(group) k <- as.numeric(group_size[as.character(group)]) const <- log_lik_pois(size, size) - log_lik - const / k + log_lik <- log_lik - const / k + log_lik[multinom_row_na(size, prob, group)] <- NA_real_ + log_lik } #' Negative Binomial Log-Likelihood diff --git a/R/params.R b/R/params.R index 502cb488..1756d916 100644 --- a/R/params.R +++ b/R/params.R @@ -18,9 +18,16 @@ #' should be returned as negative values. #' @param group A vector identifying which rows belong to the same #' multinomial trial, i.e., whose `x` values must sum to `size` and whose -#' `prob` values must sum to 1. Each group must contain at least 2 rows (a -#' trial needs at least 2 categories). Required because the categories -#' within a single multinomial trial are not independent. +#' `prob` values must sum to 1. This is for ordinary multinomial logistic +#' regression, where every trial has the same fixed set of possible +#' categories: every group must contain at least 2 rows (a trial needs at +#' least 2 categories), and every group must have the same number of rows +#' (the most common number of rows across the data). `group` must not +#' contain `NA`, since there's no way to know which trial an unlabelled +#' row belongs to. `NA` in `size` or `prob` for any row of a trial makes +#' the result `NA` for every row of that trial (not just the row it +#' appears in), since the categories within a single multinomial trial +#' are not independent -- they're drawn/scored jointly. #' @param lambda A non-negative numeric vector of means. #' @param level A number > 0 and <= 1 specifying the probability coverage of the #' interval. diff --git a/R/ran.R b/R/ran.R index dbd963c5..b35bfe19 100644 --- a/R/ran.R +++ b/R/ran.R @@ -120,7 +120,10 @@ ran_lnorm <- function(n = 1, meanlog = 0, sdlog = 1) { #' (and therefore random samples) are in \emph{long} format: one value per #' category per trial, with `group` identifying which rows belong to the #' same trial. All rows sharing a `group` must have the same `size`, and -#' their `prob` values must sum to 1. +#' their `prob` values must sum to 1. This is intended for use with +#' ordinary multinomial logistic regression, where every trial has the +#' same fixed set of possible categories: every group must have the same +#' number of rows (the most common number of rows across the data). #' #' Unlike the other `ran_*()` functions, `ran_multinom()` has no `n` #' argument: the number of samples is fully determined by `length(prob)` @@ -145,9 +148,14 @@ ran_multinom <- function(size = 1, prob, group) { size <- rep_len(size, n) prob <- rep_len(prob, n) group <- rep_len(group, n) + chk_not_any_na(group) chk_multinom_group(size, prob, group) + row_na <- multinom_row_na(size, prob, group) x <- rep(NA_real_, n) for (idx in split(seq_len(n), group)) { + if (row_na[idx[1]]) { + next + } x[idx] <- stats::rmultinom(1, size = size[idx[1]], prob = prob[idx])[, 1] } as.integer(x) diff --git a/_pkgdown.yml b/_pkgdown.yml index df964f30..630b6e77 100644 --- a/_pkgdown.yml +++ b/_pkgdown.yml @@ -106,6 +106,7 @@ reference: - '`dev_gamma_pois`' - '`dev_gamma_pois_zi`' - '`dev_lnorm`' + - '`dev_multinom`' - '`dev_neg_binom`' - '`dev_norm`' - '`dev_pois`' @@ -123,6 +124,7 @@ reference: - '`res_gamma_pois`' - '`res_gamma_pois_zi`' - '`res_lnorm`' + - '`res_multinom`' - '`res_neg_binom`' - '`res_norm`' - '`res_pois`' @@ -142,6 +144,7 @@ reference: - '`log_lik_gamma_pois_zi`' - '`log_lik_exp`' - '`log_lik_lnorm`' + - '`log_lik_multinom`' - '`log_lik_neg_binom`' - '`log_lik_norm`' - '`log_lik_pois`' @@ -199,6 +202,7 @@ reference: - '`ran_gamma_pois`' - '`ran_gamma_pois_zi`' - '`ran_lnorm`' + - '`ran_multinom`' - '`ran_neg_binom`' - '`ran_norm`' - '`ran_pois`' diff --git a/man/log_lik_multinom.Rd b/man/log_lik_multinom.Rd index 7367d3ef..7e0c485e 100644 --- a/man/log_lik_multinom.Rd +++ b/man/log_lik_multinom.Rd @@ -16,9 +16,16 @@ to 1 across the rows sharing the same \code{group}.} \item{group}{A vector identifying which rows belong to the same multinomial trial, i.e., whose \code{x} values must sum to \code{size} and whose -\code{prob} values must sum to 1. Each group must contain at least 2 rows (a -trial needs at least 2 categories). Required because the categories -within a single multinomial trial are not independent.} +\code{prob} values must sum to 1. This is for ordinary multinomial logistic +regression, where every trial has the same fixed set of possible +categories: every group must contain at least 2 rows (a trial needs at +least 2 categories), and every group must have the same number of rows +(the most common number of rows across the data). \code{group} must not +contain \code{NA}, since there's no way to know which trial an unlabelled +row belongs to. \code{NA} in \code{size} or \code{prob} for any row of a trial makes +the result \code{NA} for every row of that trial (not just the row it +appears in), since the categories within a single multinomial trial +are not independent -- they're drawn/scored jointly.} } \value{ An numeric vector of the corresponding log-likelihoods, one value diff --git a/man/params.Rd b/man/params.Rd index b2f7d92f..aad1b77f 100644 --- a/man/params.Rd +++ b/man/params.Rd @@ -21,9 +21,16 @@ should be returned as negative values.} \item{group}{A vector identifying which rows belong to the same multinomial trial, i.e., whose \code{x} values must sum to \code{size} and whose -\code{prob} values must sum to 1. Each group must contain at least 2 rows (a -trial needs at least 2 categories). Required because the categories -within a single multinomial trial are not independent.} +\code{prob} values must sum to 1. This is for ordinary multinomial logistic +regression, where every trial has the same fixed set of possible +categories: every group must contain at least 2 rows (a trial needs at +least 2 categories), and every group must have the same number of rows +(the most common number of rows across the data). \code{group} must not +contain \code{NA}, since there's no way to know which trial an unlabelled +row belongs to. \code{NA} in \code{size} or \code{prob} for any row of a trial makes +the result \code{NA} for every row of that trial (not just the row it +appears in), since the categories within a single multinomial trial +are not independent -- they're drawn/scored jointly.} \item{lambda}{A non-negative numeric vector of means.} diff --git a/man/ran_multinom.Rd b/man/ran_multinom.Rd index dd0842e7..d3a0284e 100644 --- a/man/ran_multinom.Rd +++ b/man/ran_multinom.Rd @@ -14,9 +14,16 @@ to 1 across the rows sharing the same \code{group}.} \item{group}{A vector identifying which rows belong to the same multinomial trial, i.e., whose \code{x} values must sum to \code{size} and whose -\code{prob} values must sum to 1. Each group must contain at least 2 rows (a -trial needs at least 2 categories). Required because the categories -within a single multinomial trial are not independent.} +\code{prob} values must sum to 1. This is for ordinary multinomial logistic +regression, where every trial has the same fixed set of possible +categories: every group must contain at least 2 rows (a trial needs at +least 2 categories), and every group must have the same number of rows +(the most common number of rows across the data). \code{group} must not +contain \code{NA}, since there's no way to know which trial an unlabelled +row belongs to. \code{NA} in \code{size} or \code{prob} for any row of a trial makes +the result \code{NA} for every row of that trial (not just the row it +appears in), since the categories within a single multinomial trial +are not independent -- they're drawn/scored jointly.} } \value{ An integer vector of the random samples, one per row of \code{prob}. @@ -27,7 +34,10 @@ mutually exclusive categories arising from a fixed number of trials. Data (and therefore random samples) are in \emph{long} format: one value per category per trial, with \code{group} identifying which rows belong to the same trial. All rows sharing a \code{group} must have the same \code{size}, and -their \code{prob} values must sum to 1. +their \code{prob} values must sum to 1. This is intended for use with +ordinary multinomial logistic regression, where every trial has the +same fixed set of possible categories: every group must have the same +number of rows (the most common number of rows across the data). } \details{ Unlike the other \verb{ran_*()} functions, \code{ran_multinom()} has no \code{n} diff --git a/man/res_multinom.Rd b/man/res_multinom.Rd index 15672878..9943dbf7 100644 --- a/man/res_multinom.Rd +++ b/man/res_multinom.Rd @@ -16,9 +16,16 @@ to 1 across the rows sharing the same \code{group}.} \item{group}{A vector identifying which rows belong to the same multinomial trial, i.e., whose \code{x} values must sum to \code{size} and whose -\code{prob} values must sum to 1. Each group must contain at least 2 rows (a -trial needs at least 2 categories). Required because the categories -within a single multinomial trial are not independent.} +\code{prob} values must sum to 1. This is for ordinary multinomial logistic +regression, where every trial has the same fixed set of possible +categories: every group must contain at least 2 rows (a trial needs at +least 2 categories), and every group must have the same number of rows +(the most common number of rows across the data). \code{group} must not +contain \code{NA}, since there's no way to know which trial an unlabelled +row belongs to. \code{NA} in \code{size} or \code{prob} for any row of a trial makes +the result \code{NA} for every row of that trial (not just the row it +appears in), since the categories within a single multinomial trial +are not independent -- they're drawn/scored jointly.} \item{type}{A string of the residual type. 'raw' for raw residuals 'dev' for deviance residuals and 'data' for the data.} diff --git a/tests/testthat/test-log-lik.R b/tests/testthat/test-log-lik.R index 01549db9..90fb0038 100644 --- a/tests/testthat/test-log-lik.R +++ b/tests/testthat/test-log-lik.R @@ -349,6 +349,49 @@ test_that("log_lik_multinom", { log_lik_multinom(4, 10, 1, 1), "must contain at least 2 rows" ) + # group must not contain NA -- there's no way to know which trial an + # unlabelled row belongs to + expect_error( + log_lik_multinom(c(4, 6, 10), c(10, 10, 10), c(0.4, 0.6, 1), c(1, 1, NA)), + "must not have any missing values" + ) + expect_error( + log_lik_multinom(c(4, 6), c(10, 10), c(0.4, 0.6), c(NA, NA)), + "must not have any missing values" + ) + # NA in size or prob is a structural input shared by the whole trial, so + # it makes the whole trial's result NA, not just the row where it appears + expect_identical( + log_lik_multinom(c(4, 6), c(10, NA), c(0.4, 0.6), c(1, 1)), + c(NA_real_, NA_real_) + ) + # NA in x is just a missing observation for that one category -- it + # doesn't need the rest of the trial to compute its own log_lik, so it + # doesn't taint sibling rows the way NA size/prob/group does + res <- log_lik_multinom(c(4, NA), c(10, 10), c(0.4, 0.6), c(1, 1)) + expect_false(is.na(res[1])) + expect_true(is.na(res[2])) + # an NA elsewhere doesn't leak into an unrelated, fully-known group (same + # number of categories in both groups, since groups must match on that) + res <- log_lik_multinom( + c(4, 3, 3, 2, 2, 2), + c(10, 10, NA, 6, 6, 6), + c(0.2, 0.3, 0.5, 0.2, 0.3, 0.5), + c(1, 1, 1, 2, 2, 2) + ) + expect_identical(res[1:3], c(NA_real_, NA_real_, NA_real_)) + expect_false(anyNA(res[4:6])) + # every group must have the same number of rows (categories) as the most + # common number of rows per group in the data + expect_error( + log_lik_multinom( + c(4, 3, 3, 5, 5, 3, 3), + c(10, 10, 10, 10, 10, 6, 6), + c(0.2, 0.3, 0.5, 0.5, 0.5, 0.5, 0.5), + c(1, 1, 1, 2, 2, 3, 3) + ), + "Every `group` should have the same number of rows" + ) expect_equal( sum(log_lik_multinom( c(1, 3, 6), @@ -358,11 +401,11 @@ test_that("log_lik_multinom", { )), dmultinom(c(1, 3, 6), size = 10, prob = c(0.2, 0.3, 0.5), log = TRUE) ) - # multiple trials in long format - x <- c(1, 3, 6, 2, 2) - size <- c(10, 10, 10, 4, 4) - prob <- c(0.2, 0.3, 0.5, 0.5, 0.5) - group <- c(1, 1, 1, 2, 2) + # multiple trials in long format (same number of categories per group) + x <- c(1, 3, 6, 2, 1, 1) + size <- c(10, 10, 10, 4, 4, 4) + prob <- c(0.2, 0.3, 0.5, 0.5, 0.25, 0.25) + group <- c(1, 1, 1, 2, 2, 2) ll <- log_lik_multinom(x, size, prob, group) expect_equal( sum(ll[group == 1]), diff --git a/tests/testthat/test-ran.R b/tests/testthat/test-ran.R index be711fcf..12afa41c 100644 --- a/tests/testthat/test-ran.R +++ b/tests/testthat/test-ran.R @@ -132,6 +132,49 @@ test_that("ran_multinom", { ), "must contain at least 2 rows" ) + # every group must have the same number of rows (categories) as the most + # common number of rows per group in the data (ordinary multinomial + # logistic regression assumes a fixed set of categories for every trial) + expect_error( + ran_multinom( + size = c(10, 10, 10, 10, 10, 10, 6, 6), + prob = c(0.2, 0.3, 0.5, 0.2, 0.3, 0.5, 0.5, 0.5), + group = c(1, 1, 1, 2, 2, 2, 3, 3) + ), + "Every `group` should have the same number of rows" + ) + # group must not contain NA -- there's no way to know which trial an + # unlabelled row belongs to + expect_error( + ran_multinom(size = c(10, 10, 10), prob = c(0.4, 0.6, 1), group = c(1, 1, NA)), + "must not have any missing values" + ) + expect_error( + ran_multinom(size = c(10, 10), prob = c(0.4, 0.6), group = c(NA, NA)), + "must not have any missing values" + ) + # NA in size or prob is a structural input shared by the whole trial (the + # categories aren't independent draws), so it makes the whole trial's + # result NA, not just the row where the NA appears + expect_identical( + ran_multinom(size = c(10, NA), prob = c(0.4, 0.6), group = c(1, 1)), + c(NA_integer_, NA_integer_) + ) + expect_identical( + ran_multinom(size = c(10, 10), prob = c(0.4, NA), group = c(1, 1)), + c(NA_integer_, NA_integer_) + ) + # an NA elsewhere doesn't leak into an unrelated, fully-known group (same + # number of categories in both groups, since groups must match on that) + withr::with_seed(101, { + x <- ran_multinom( + size = c(10, 10, NA, 6, 6, 6), + prob = c(0.2, 0.3, 0.5, 0.2, 0.3, 0.5), + group = c(1, 1, 1, 2, 2, 2) + ) + expect_identical(x[1:3], c(NA_integer_, NA_integer_, NA_integer_)) + expect_identical(sum(x[4:6]), 6L) + }) withr::with_seed(101, { x <- ran_multinom(size = 10, prob = c(0.2, 0.3, 0.5), group = c(1, 1, 1)) expect_identical(sum(x), 10L) diff --git a/tests/testthat/test-res.R b/tests/testthat/test-res.R index 715f66f5..3d615875 100644 --- a/tests/testthat/test-res.R +++ b/tests/testthat/test-res.R @@ -536,6 +536,49 @@ test_that("res_multinom", { res_multinom(4, 10, 1, 1, simulate = TRUE), "must contain at least 2 rows" ) + # NA in size or prob is a structural input shared by the whole trial, so + # it makes the whole trial's result NA, not just the row where it appears + expect_identical( + res_multinom(c(4, 6), c(10, NA), c(0.4, 0.6), c(1, 1), type = "data", simulate = TRUE), + c(NA_integer_, NA_integer_) + ) + # every group must have the same number of rows (categories) as the most + # common number of rows per group in the data + expect_error( + res_multinom( + c(4, 3, 3, 5, 5), + c(10, 10, 10, 6, 6), + c(0.2, 0.3, 0.5, 0.5, 0.5), + c(1, 1, 1, 2, 2), + type = "data", + simulate = TRUE + ), + "Every `group` should have the same number of rows" + ) + # group must not contain NA -- there's no way to know which trial an + # unlabelled row belongs to + expect_error( + res_multinom( + c(4, 6, 10), + c(10, 10, 10), + c(0.4, 0.6, 1), + c(1, 1, NA), + type = "data", + simulate = TRUE + ), + "must not have any missing values" + ) + expect_error( + res_multinom( + c(4, 6), + c(10, 10), + c(0.4, 0.6), + c(NA, NA), + type = "data", + simulate = TRUE + ), + "must not have any missing values" + ) # sum of squared deviance residuals recovers the row-level deviance expect_equal( From 3bd467dd92ea6c64ead7e1d59d79b5b9f5b405c2 Mon Sep 17 00:00:00 2001 From: Nicole Hill Date: Fri, 7 Aug 2026 11:01:16 -0700 Subject: [PATCH 6/9] Condense multinomial docs/comments; validate group in dev_multinom and res_multinom(simulate = FALSE) Trim verbose internal comments and roxygen prose added for the multinomial functions, and move NA-propagation notes to the size/prob docs where they actually apply. dev_multinom() now takes group and validates size/prob against it, and res_multinom() validates group even when simulate = FALSE, so a prob vector that doesn't sum to 1 can no longer silently produce a wrong result on either path. --- R/dev.R | 38 ++++++++++++++--------------- R/internal.R | 50 ++++++++------------------------------- R/log-lik.R | 31 ++++++++++++------------ R/params.R | 15 ++++-------- R/ran.R | 20 +++++++--------- R/res.R | 39 ++++++++++++++++-------------- man/dev_multinom.Rd | 36 ++++++++++++++-------------- man/log_lik_multinom.Rd | 44 ++++++++++++++-------------------- man/params.Rd | 15 ++++-------- man/ran_multinom.Rd | 35 ++++++++++----------------- man/res_multinom.Rd | 44 +++++++++++++--------------------- tests/testthat/test-dev.R | 30 ++++++++++++++++------- tests/testthat/test-res.R | 12 +++++++++- 13 files changed, 181 insertions(+), 228 deletions(-) diff --git a/R/dev.R b/R/dev.R index 7bca89e7..38b2380e 100644 --- a/R/dev.R +++ b/R/dev.R @@ -269,35 +269,35 @@ dev_lnorm <- function(x, meanlog = 0, sdlog = 1, res = FALSE) { #' Multinomial Deviances #' -#' The multinomial distribution models the counts across two or more -#' mutually exclusive categories arising from a fixed number of trials. Data -#' are in \emph{long} format: one row per category per trial (see -#' [log_lik_multinom()] for details of the `group` argument used elsewhere -#' in the multinomial family). -#' -#' Unlike the other `dev_*()` functions, `dev_multinom()` doesn't take a -#' `group` argument, because a category's deviance contribution only depends -#' on its own `x` and `mu = size * prob` -- the multinomial coefficient -#' cancels out of the deviance difference. This is the Poisson-equivalent -#' deviance (see [dev_pois()]): summing it over the rows of one trial -#' recovers the trial's exact multinomial deviance, provided `prob` sums to -#' 1 across those rows. Because each cell's deviance is an ordinary scalar -#' comparison of `x` to `mu`, its deviance residual has an ordinary sign -- -#' there's no need for the "extended sign" methods proposed elsewhere for -#' multinomial residuals. +#' Models the counts across two or more mutually exclusive categories from a +#' fixed number of trials, in \emph{long} format: one row per category per +#' trial, with `group` identifying which rows belong to the same trial. +#' +#' A category's deviance depends only on its own `x` and `mu = size * prob`, +#' not on the rest of its trial, so `group` is used only to validate `size` +#' and `prob` (see [log_lik_multinom()]), not in the calculation itself. +#' `dev_multinom()` is the Poisson-equivalent deviance (see [dev_pois()]): +#' summing it over a trial's rows recovers the trial's exact multinomial +#' deviance. #' #' @inheritParams params #' @param x A non-negative whole numeric vector of the category counts. #' @param prob A numeric vector of the probability of the category. Must sum -#' to 1 across the rows belonging to the same trial. +#' to 1 across the rows sharing the same `group`. #' #' @return An numeric vector of the corresponding deviances or deviance residuals. #' @family dev_dist #' @export #' #' @examples -#' dev_multinom(c(1, 3, 6), size = 10, prob = c(0.2, 0.3, 0.5)) -dev_multinom <- function(x, size = 1, prob, res = FALSE) { +#' dev_multinom(c(1, 3, 6), size = 10, prob = c(0.2, 0.3, 0.5), group = c(1, 1, 1)) +dev_multinom <- function(x, size = 1, prob, group, res = FALSE) { + n <- length(x) + size <- rep_len(size, n) + prob <- rep_len(prob, n) + group <- rep_len(group, n) + chk_not_any_na(group) + chk_multinom_group(size, prob, group) dev_pois(x, lambda = size * prob, res = res) } diff --git a/R/internal.R b/R/internal.R index 9d3bac5c..94788643 100644 --- a/R/internal.R +++ b/R/internal.R @@ -2,39 +2,14 @@ dev_res <- function(x, mu, dev) { sign(x - mu) * sqrt(dev) } -# Validates that every row sharing a `group` (multinomial trial) has the -# same `size` and that its `prob` values sum to 1, since both the -# Poisson-equivalent deviance/log-likelihood identities and simulation via -# rmultinom() require this to hold. Also rejects singleton groups: a group -# with only one row can't represent a real multinomial trial (a trial needs -# at least 2 categories) -- rmultinom() would just deterministically return -# `size` for it. This typically happens if the function ends up being -# evaluated separately for each row instead of once for the whole data -# vector, so each call only ever sees a single row's own group value. -# -# Callers must reject NA `group` values (chk_not_any_na(group)) before -# calling this: unlike `size`/`prob` (attributes of an already-identified -# trial, where a missing value leaves that one trial undefined -- see -# multinom_row_na()), `group` is what identifies the trial in the first -# place. There's no sensible way to know which trial (if any) an unlabelled -# row belongs to, so it can't be resolved by attributing NA to "the trial it -# belongs to" the way size/prob can. That check lives in each calling -# function rather than here so it runs before `group` is used to split the -# data below. -# -# The size/prob checks compare only the non-NA values within a group, so a -# merely-missing value doesn't trigger a spurious error -- it's handled by -# multinom_row_na() instead, which marks the whole group's output NA rather -# than erroring (see its own comment for why). -# -# Also compares every group's row count (its number of categories) against -# the most common row count across all groups in the call, erroring on any -# that differ. This is specifically for ordinary multinomial logistic -# regression, where every trial has the same fixed set of possible -# categories: a group with fewer rows than the rest of the dataset is -# almost always a trial silently missing a row, which would otherwise be -# validated/simulated/scored as if it were a genuinely smaller, complete -# trial, giving a wrong (not NA, not an error) answer. +# Checks every group shares one `size` and `prob` values summing to 1 +# (required by rmultinom() and the deviance/log-lik identities), has >= 2 +# rows (a trial needs >= 2 categories -- singletons usually mean `group` +# was evaluated row-by-row instead of over the whole vector), and matches +# the modal row count across groups (a short group usually means a row was +# lost). Only non-NA values are compared, so lone NAs don't error here -- +# see multinom_row_na(). Callers must chk_not_any_na(group) first; `group` +# itself can't be NA-tolerant since it's what identifies the trial. chk_multinom_group <- function(size, prob, group) { groups <- split(seq_along(group), group) for (idx in groups) { @@ -78,13 +53,8 @@ chk_multinom_group <- function(size, prob, group) { } } -# Returns, for each row, whether its multinomial trial's draw/score can't be -# computed because `size`/`prob` is NA for ANY row sharing its group. A -# multinomial trial is one joint computation across all of its rows (the -# categories aren't independent), so missing information anywhere in that -# trial leaves the whole trial's result undefined -- not just the row that -# happens to be missing a value. `group` itself can't be NA (rejected by -# chk_multinom_group()), so every row here belongs to a real, known trial. +# Flags every row whose trial has an NA `size`/`prob` anywhere in the group, +# since a trial's categories are scored/drawn jointly, not independently. multinom_row_na <- function(size, prob, group) { bad <- is.na(size) | is.na(prob) result <- rep(FALSE, length(group)) diff --git a/R/log-lik.R b/R/log-lik.R index f5c646d5..cf88645e 100644 --- a/R/log-lik.R +++ b/R/log-lik.R @@ -429,26 +429,25 @@ log_lik_lnorm <- function(x, meanlog = 0, sdlog = 1, tlower = 0, tupper = Inf) { #' Multinomial Log-Likelihood #' -#' The multinomial distribution models the counts across two or more -#' mutually exclusive categories arising from a fixed number of trials. Data -#' are in \emph{long} format: one row per category per trial, with `group` -#' identifying which rows belong to the same trial. All rows sharing a -#' `group` must have the same `size`, and their `prob` values must sum to 1. -#' -#' A trial's log-likelihood doesn't split evenly across its category rows, -#' because the multinomial coefficient is a property of the whole trial, not -#' any one category. `log_lik_multinom()` uses the identity that a -#' multinomial is equivalent to independent Poissons conditional on the -#' trial total: the log-likelihood of category `k` of trial `i` is the -#' Poisson log-likelihood of `x` given `mu = size * prob`, minus an even -#' share of the trial's normalizing constant (so that summing -#' `log_lik_multinom()` over the rows of one `group` recovers the trial's -#' exact multinomial log-likelihood). +#' Models the counts across two or more mutually exclusive categories from a +#' fixed number of trials, in \emph{long} format: one row per category per +#' trial, with `group` identifying which rows belong to the same trial. All +#' rows sharing a `group` must have the same `size`, and their `prob` values +#' must sum to 1. +#' +#' A trial's log-likelihood doesn't split evenly across its rows, since the +#' multinomial coefficient belongs to the whole trial. `log_lik_multinom()` +#' uses the multinomial-as-independent-Poissons identity: each row's value +#' is the Poisson log-likelihood of `x` given `mu = size * prob`, minus an +#' even share of the trial's normalizing constant, so summing over a +#' `group` recovers the trial's exact multinomial log-likelihood. #' #' @inheritParams params #' @param x A non-negative whole numeric vector of the category counts. #' @param prob A numeric vector of the probability of the category. Must sum -#' to 1 across the rows sharing the same `group`. +#' to 1 across the rows sharing the same `group`. `NA` in `size` or `prob` +#' for any row of a trial makes the log-likelihood `NA` for every row of +#' that trial, since a trial's categories are scored jointly. #' #' @return An numeric vector of the corresponding log-likelihoods, one value #' per row of `x`. diff --git a/R/params.R b/R/params.R index 1756d916..357c8ab5 100644 --- a/R/params.R +++ b/R/params.R @@ -17,17 +17,10 @@ #' @param directional A flag specifying whether probabilities less than 0.5 #' should be returned as negative values. #' @param group A vector identifying which rows belong to the same -#' multinomial trial, i.e., whose `x` values must sum to `size` and whose -#' `prob` values must sum to 1. This is for ordinary multinomial logistic -#' regression, where every trial has the same fixed set of possible -#' categories: every group must contain at least 2 rows (a trial needs at -#' least 2 categories), and every group must have the same number of rows -#' (the most common number of rows across the data). `group` must not -#' contain `NA`, since there's no way to know which trial an unlabelled -#' row belongs to. `NA` in `size` or `prob` for any row of a trial makes -#' the result `NA` for every row of that trial (not just the row it -#' appears in), since the categories within a single multinomial trial -#' are not independent -- they're drawn/scored jointly. +#' multinomial trial (whose `x` values sum to `size` and `prob` values sum +#' to 1). Every group must have at least 2 rows and the same number of +#' rows as the rest of the data (a fixed set of categories, as in +#' multinomial logistic regression), and must not contain `NA`. #' @param lambda A non-negative numeric vector of means. #' @param level A number > 0 and <= 1 specifying the probability coverage of the #' interval. diff --git a/R/ran.R b/R/ran.R index b35bfe19..e6bb6a8f 100644 --- a/R/ran.R +++ b/R/ran.R @@ -115,24 +115,22 @@ ran_lnorm <- function(n = 1, meanlog = 0, sdlog = 1) { #' Multinomial Random Samples #' -#' The multinomial distribution models the counts across two or more -#' mutually exclusive categories arising from a fixed number of trials. Data -#' (and therefore random samples) are in \emph{long} format: one value per -#' category per trial, with `group` identifying which rows belong to the -#' same trial. All rows sharing a `group` must have the same `size`, and -#' their `prob` values must sum to 1. This is intended for use with -#' ordinary multinomial logistic regression, where every trial has the -#' same fixed set of possible categories: every group must have the same -#' number of rows (the most common number of rows across the data). +#' Models the counts across two or more mutually exclusive categories from a +#' fixed number of trials, in \emph{long} format: one value per category per +#' trial, with `group` identifying which rows belong to the same trial. All +#' rows sharing a `group` must have the same `size`, and their `prob` values +#' must sum to 1. #' #' Unlike the other `ran_*()` functions, `ran_multinom()` has no `n` #' argument: the number of samples is fully determined by `length(prob)` -#' (equivalently `length(group)`), because a trial's categories can't be +#' (equivalently `length(group)`), since a trial's categories can't be #' generated independently of one another. #' #' @inheritParams params #' @param prob A numeric vector of the probability of the category. Must sum -#' to 1 across the rows sharing the same `group`. +#' to 1 across the rows sharing the same `group`. `NA` in `size` or `prob` +#' for any row of a trial makes the sample `NA` for every row of that +#' trial, since a trial's categories are drawn jointly. #' @return An integer vector of the random samples, one per row of `prob`. #' @family ran_dist #' @export diff --git a/R/res.R b/R/res.R index 02a46df4..a4f4f4a4 100644 --- a/R/res.R +++ b/R/res.R @@ -253,24 +253,21 @@ res_lnorm <- function( #' Multinomial Residuals #' -#' The multinomial distribution models the counts across two or more -#' mutually exclusive categories arising from a fixed number of trials. Data -#' are in \emph{long} format: one row per category per trial, with `group` -#' identifying which rows belong to the same trial (see +#' Models the counts across two or more mutually exclusive categories from a +#' fixed number of trials, in \emph{long} format: one row per category per +#' trial, with `group` identifying which rows belong to the same trial (see #' [log_lik_multinom()] for details). `res_multinom()` returns one residual -#' per row (per category per trial), not one per trial -- since a trial's -#' categories are not independent, there's no single meaningful residual (or -#' sign) for the trial as a whole. The classic per-trial deviance statistic -#' can always be recovered by summing the squared `type = "dev"` residuals -#' within a `group`. -#' -#' `group` is only used when `simulate = TRUE`, to draw a joint, -#' correlation-preserving replicate for each trial (via [ran_multinom()]) -#' rather than simulating each category independently, which would -#' understate the true (negative) covariance among a trial's categories. -#' This requires `res_multinom()` to be evaluated on the full data vector at -#' once, so that every row sharing a `group` is visible in the same call, -#' rather than evaluated separately for each row. +#' per row, not one per trial, since a trial's categories aren't +#' independent and so have no single meaningful residual as a whole; the +#' classic per-trial deviance statistic can be recovered by summing the +#' squared `type = "dev"` residuals within a `group`. +#' +#' `group` is validated (same `size`, `prob` summing to 1, no singleton or +#' short groups, no `NA`) regardless of `simulate`, but is only otherwise +#' used when `simulate = TRUE`, to draw a joint, correlation-preserving +#' replicate per trial (via [ran_multinom()]) rather than simulating each +#' category independently, which requires `res_multinom()` to see every row +#' of a `group` in the same call. #' #' @inheritParams params #' @param x A non-negative whole numeric vector of the category counts. @@ -292,6 +289,12 @@ res_multinom <- function( simulate = FALSE ) { chk_string(type) + n <- length(x) + size <- rep_len(size, n) + prob <- rep_len(prob, n) + group <- rep_len(group, n) + chk_not_any_na(group) + chk_multinom_group(size, prob, group) if (!vld_false(simulate)) { x <- ran_multinom(size = size, prob = prob, group = group) } @@ -301,7 +304,7 @@ res_multinom <- function( data = x, raw = x - mu, standardized = (x - mu) / sqrt(mu * (1 - prob)), - dev = dev_multinom(x, size = size, prob = prob, res = TRUE), + dev = dev_multinom(x, size = size, prob = prob, group = group, res = TRUE), chk_subset(x, c("data", "raw", "dev", "standardized")) ) } diff --git a/man/dev_multinom.Rd b/man/dev_multinom.Rd index 472ac5e1..57dd2172 100644 --- a/man/dev_multinom.Rd +++ b/man/dev_multinom.Rd @@ -4,7 +4,7 @@ \alias{dev_multinom} \title{Multinomial Deviances} \usage{ -dev_multinom(x, size = 1, prob, res = FALSE) +dev_multinom(x, size = 1, prob, group, res = FALSE) } \arguments{ \item{x}{A non-negative whole numeric vector of the category counts.} @@ -12,7 +12,13 @@ dev_multinom(x, size = 1, prob, res = FALSE) \item{size}{A non-negative whole numeric vector of the number of trials.} \item{prob}{A numeric vector of the probability of the category. Must sum -to 1 across the rows belonging to the same trial.} +to 1 across the rows sharing the same \code{group}.} + +\item{group}{A vector identifying which rows belong to the same +multinomial trial (whose \code{x} values sum to \code{size} and \code{prob} values sum +to 1). Every group must have at least 2 rows and the same number of +rows as the rest of the data (a fixed set of categories, as in +multinomial logistic regression), and must not contain \code{NA}.} \item{res}{A flag specifying whether to return the deviance residual as opposed to the deviance.} @@ -21,26 +27,20 @@ opposed to the deviance.} An numeric vector of the corresponding deviances or deviance residuals. } \description{ -The multinomial distribution models the counts across two or more -mutually exclusive categories arising from a fixed number of trials. Data -are in \emph{long} format: one row per category per trial (see -\code{\link[=log_lik_multinom]{log_lik_multinom()}} for details of the \code{group} argument used elsewhere -in the multinomial family). +Models the counts across two or more mutually exclusive categories from a +fixed number of trials, in \emph{long} format: one row per category per +trial, with \code{group} identifying which rows belong to the same trial. } \details{ -Unlike the other \verb{dev_*()} functions, \code{dev_multinom()} doesn't take a -\code{group} argument, because a category's deviance contribution only depends -on its own \code{x} and \code{mu = size * prob} -- the multinomial coefficient -cancels out of the deviance difference. This is the Poisson-equivalent -deviance (see \code{\link[=dev_pois]{dev_pois()}}): summing it over the rows of one trial -recovers the trial's exact multinomial deviance, provided \code{prob} sums to -1 across those rows. Because each cell's deviance is an ordinary scalar -comparison of \code{x} to \code{mu}, its deviance residual has an ordinary sign -- -there's no need for the "extended sign" methods proposed elsewhere for -multinomial residuals. +A category's deviance depends only on its own \code{x} and \code{mu = size * prob}, +not on the rest of its trial, so \code{group} is used only to validate \code{size} +and \code{prob} (see \code{\link[=log_lik_multinom]{log_lik_multinom()}}), not in the calculation itself. +\code{dev_multinom()} is the Poisson-equivalent deviance (see \code{\link[=dev_pois]{dev_pois()}}): +summing it over a trial's rows recovers the trial's exact multinomial +deviance. } \examples{ -dev_multinom(c(1, 3, 6), size = 10, prob = c(0.2, 0.3, 0.5)) +dev_multinom(c(1, 3, 6), size = 10, prob = c(0.2, 0.3, 0.5), group = c(1, 1, 1)) } \seealso{ Other dev_dist: diff --git a/man/log_lik_multinom.Rd b/man/log_lik_multinom.Rd index 7e0c485e..f7b1d8f5 100644 --- a/man/log_lik_multinom.Rd +++ b/man/log_lik_multinom.Rd @@ -12,42 +12,34 @@ log_lik_multinom(x, size = 1, prob, group) \item{size}{A non-negative whole numeric vector of the number of trials.} \item{prob}{A numeric vector of the probability of the category. Must sum -to 1 across the rows sharing the same \code{group}.} +to 1 across the rows sharing the same \code{group}. \code{NA} in \code{size} or \code{prob} +for any row of a trial makes the log-likelihood \code{NA} for every row of +that trial, since a trial's categories are scored jointly.} \item{group}{A vector identifying which rows belong to the same -multinomial trial, i.e., whose \code{x} values must sum to \code{size} and whose -\code{prob} values must sum to 1. This is for ordinary multinomial logistic -regression, where every trial has the same fixed set of possible -categories: every group must contain at least 2 rows (a trial needs at -least 2 categories), and every group must have the same number of rows -(the most common number of rows across the data). \code{group} must not -contain \code{NA}, since there's no way to know which trial an unlabelled -row belongs to. \code{NA} in \code{size} or \code{prob} for any row of a trial makes -the result \code{NA} for every row of that trial (not just the row it -appears in), since the categories within a single multinomial trial -are not independent -- they're drawn/scored jointly.} +multinomial trial (whose \code{x} values sum to \code{size} and \code{prob} values sum +to 1). Every group must have at least 2 rows and the same number of +rows as the rest of the data (a fixed set of categories, as in +multinomial logistic regression), and must not contain \code{NA}.} } \value{ An numeric vector of the corresponding log-likelihoods, one value per row of \code{x}. } \description{ -The multinomial distribution models the counts across two or more -mutually exclusive categories arising from a fixed number of trials. Data -are in \emph{long} format: one row per category per trial, with \code{group} -identifying which rows belong to the same trial. All rows sharing a -\code{group} must have the same \code{size}, and their \code{prob} values must sum to 1. +Models the counts across two or more mutually exclusive categories from a +fixed number of trials, in \emph{long} format: one row per category per +trial, with \code{group} identifying which rows belong to the same trial. All +rows sharing a \code{group} must have the same \code{size}, and their \code{prob} values +must sum to 1. } \details{ -A trial's log-likelihood doesn't split evenly across its category rows, -because the multinomial coefficient is a property of the whole trial, not -any one category. \code{log_lik_multinom()} uses the identity that a -multinomial is equivalent to independent Poissons conditional on the -trial total: the log-likelihood of category \code{k} of trial \code{i} is the -Poisson log-likelihood of \code{x} given \code{mu = size * prob}, minus an even -share of the trial's normalizing constant (so that summing -\code{log_lik_multinom()} over the rows of one \code{group} recovers the trial's -exact multinomial log-likelihood). +A trial's log-likelihood doesn't split evenly across its rows, since the +multinomial coefficient belongs to the whole trial. \code{log_lik_multinom()} +uses the multinomial-as-independent-Poissons identity: each row's value +is the Poisson log-likelihood of \code{x} given \code{mu = size * prob}, minus an +even share of the trial's normalizing constant, so summing over a +\code{group} recovers the trial's exact multinomial log-likelihood. } \examples{ log_lik_multinom(c(1, 3, 6), size = 10, prob = c(0.2, 0.3, 0.5), group = c(1, 1, 1)) diff --git a/man/params.Rd b/man/params.Rd index aad1b77f..4105f9d4 100644 --- a/man/params.Rd +++ b/man/params.Rd @@ -20,17 +20,10 @@ level.} should be returned as negative values.} \item{group}{A vector identifying which rows belong to the same -multinomial trial, i.e., whose \code{x} values must sum to \code{size} and whose -\code{prob} values must sum to 1. This is for ordinary multinomial logistic -regression, where every trial has the same fixed set of possible -categories: every group must contain at least 2 rows (a trial needs at -least 2 categories), and every group must have the same number of rows -(the most common number of rows across the data). \code{group} must not -contain \code{NA}, since there's no way to know which trial an unlabelled -row belongs to. \code{NA} in \code{size} or \code{prob} for any row of a trial makes -the result \code{NA} for every row of that trial (not just the row it -appears in), since the categories within a single multinomial trial -are not independent -- they're drawn/scored jointly.} +multinomial trial (whose \code{x} values sum to \code{size} and \code{prob} values sum +to 1). Every group must have at least 2 rows and the same number of +rows as the rest of the data (a fixed set of categories, as in +multinomial logistic regression), and must not contain \code{NA}.} \item{lambda}{A non-negative numeric vector of means.} diff --git a/man/ran_multinom.Rd b/man/ran_multinom.Rd index d3a0284e..de3bcc04 100644 --- a/man/ran_multinom.Rd +++ b/man/ran_multinom.Rd @@ -10,39 +10,30 @@ ran_multinom(size = 1, prob, group) \item{size}{A non-negative whole numeric vector of the number of trials.} \item{prob}{A numeric vector of the probability of the category. Must sum -to 1 across the rows sharing the same \code{group}.} +to 1 across the rows sharing the same \code{group}. \code{NA} in \code{size} or \code{prob} +for any row of a trial makes the sample \code{NA} for every row of that +trial, since a trial's categories are drawn jointly.} \item{group}{A vector identifying which rows belong to the same -multinomial trial, i.e., whose \code{x} values must sum to \code{size} and whose -\code{prob} values must sum to 1. This is for ordinary multinomial logistic -regression, where every trial has the same fixed set of possible -categories: every group must contain at least 2 rows (a trial needs at -least 2 categories), and every group must have the same number of rows -(the most common number of rows across the data). \code{group} must not -contain \code{NA}, since there's no way to know which trial an unlabelled -row belongs to. \code{NA} in \code{size} or \code{prob} for any row of a trial makes -the result \code{NA} for every row of that trial (not just the row it -appears in), since the categories within a single multinomial trial -are not independent -- they're drawn/scored jointly.} +multinomial trial (whose \code{x} values sum to \code{size} and \code{prob} values sum +to 1). Every group must have at least 2 rows and the same number of +rows as the rest of the data (a fixed set of categories, as in +multinomial logistic regression), and must not contain \code{NA}.} } \value{ An integer vector of the random samples, one per row of \code{prob}. } \description{ -The multinomial distribution models the counts across two or more -mutually exclusive categories arising from a fixed number of trials. Data -(and therefore random samples) are in \emph{long} format: one value per -category per trial, with \code{group} identifying which rows belong to the -same trial. All rows sharing a \code{group} must have the same \code{size}, and -their \code{prob} values must sum to 1. This is intended for use with -ordinary multinomial logistic regression, where every trial has the -same fixed set of possible categories: every group must have the same -number of rows (the most common number of rows across the data). +Models the counts across two or more mutually exclusive categories from a +fixed number of trials, in \emph{long} format: one value per category per +trial, with \code{group} identifying which rows belong to the same trial. All +rows sharing a \code{group} must have the same \code{size}, and their \code{prob} values +must sum to 1. } \details{ Unlike the other \verb{ran_*()} functions, \code{ran_multinom()} has no \code{n} argument: the number of samples is fully determined by \code{length(prob)} -(equivalently \code{length(group)}), because a trial's categories can't be +(equivalently \code{length(group)}), since a trial's categories can't be generated independently of one another. } \examples{ diff --git a/man/res_multinom.Rd b/man/res_multinom.Rd index 9943dbf7..ce90ef38 100644 --- a/man/res_multinom.Rd +++ b/man/res_multinom.Rd @@ -15,17 +15,10 @@ res_multinom(x, size = 1, prob, group, type = "dev", simulate = FALSE) to 1 across the rows sharing the same \code{group}.} \item{group}{A vector identifying which rows belong to the same -multinomial trial, i.e., whose \code{x} values must sum to \code{size} and whose -\code{prob} values must sum to 1. This is for ordinary multinomial logistic -regression, where every trial has the same fixed set of possible -categories: every group must contain at least 2 rows (a trial needs at -least 2 categories), and every group must have the same number of rows -(the most common number of rows across the data). \code{group} must not -contain \code{NA}, since there's no way to know which trial an unlabelled -row belongs to. \code{NA} in \code{size} or \code{prob} for any row of a trial makes -the result \code{NA} for every row of that trial (not just the row it -appears in), since the categories within a single multinomial trial -are not independent -- they're drawn/scored jointly.} +multinomial trial (whose \code{x} values sum to \code{size} and \code{prob} values sum +to 1). Every group must have at least 2 rows and the same number of +rows as the rest of the data (a fixed set of categories, as in +multinomial logistic regression), and must not contain \code{NA}.} \item{type}{A string of the residual type. 'raw' for raw residuals 'dev' for deviance residuals and 'data' for the data.} @@ -36,25 +29,22 @@ deviance residuals and 'data' for the data.} An numeric vector of the corresponding residuals. } \description{ -The multinomial distribution models the counts across two or more -mutually exclusive categories arising from a fixed number of trials. Data -are in \emph{long} format: one row per category per trial, with \code{group} -identifying which rows belong to the same trial (see +Models the counts across two or more mutually exclusive categories from a +fixed number of trials, in \emph{long} format: one row per category per +trial, with \code{group} identifying which rows belong to the same trial (see \code{\link[=log_lik_multinom]{log_lik_multinom()}} for details). \code{res_multinom()} returns one residual -per row (per category per trial), not one per trial -- since a trial's -categories are not independent, there's no single meaningful residual (or -sign) for the trial as a whole. The classic per-trial deviance statistic -can always be recovered by summing the squared \code{type = "dev"} residuals -within a \code{group}. +per row, not one per trial, since a trial's categories aren't +independent and so have no single meaningful residual as a whole; the +classic per-trial deviance statistic can be recovered by summing the +squared \code{type = "dev"} residuals within a \code{group}. } \details{ -\code{group} is only used when \code{simulate = TRUE}, to draw a joint, -correlation-preserving replicate for each trial (via \code{\link[=ran_multinom]{ran_multinom()}}) -rather than simulating each category independently, which would -understate the true (negative) covariance among a trial's categories. -This requires \code{res_multinom()} to be evaluated on the full data vector at -once, so that every row sharing a \code{group} is visible in the same call, -rather than evaluated separately for each row. +\code{group} is validated (same \code{size}, \code{prob} summing to 1, no singleton or +short groups, no \code{NA}) regardless of \code{simulate}, but is only otherwise +used when \code{simulate = TRUE}, to draw a joint, correlation-preserving +replicate per trial (via \code{\link[=ran_multinom]{ran_multinom()}}) rather than simulating each +category independently, which requires \code{res_multinom()} to see every row +of a \code{group} in the same call. } \examples{ res_multinom(c(1, 3, 6), size = 10, prob = c(0.2, 0.3, 0.5), group = c(1, 1, 1)) diff --git a/tests/testthat/test-dev.R b/tests/testthat/test-dev.R index 8351d1a8..784e5a81 100644 --- a/tests/testthat/test-dev.R +++ b/tests/testthat/test-dev.R @@ -475,26 +475,37 @@ test_that("dev_lnorm deviance", { test_that("dev_multinom", { expect_identical( - dev_multinom(numeric(0), numeric(0), numeric(0)), + dev_multinom(numeric(0), numeric(0), numeric(0), group = numeric(0)), numeric(0) ) expect_identical( - dev_multinom(c(2, 3, 5), size = 10, prob = c(0.2, 0.3, 0.5)), + dev_multinom(c(2, 3, 5), size = 10, prob = c(0.2, 0.3, 0.5), group = c(1, 1, 1)), c(0, 0, 0) ) expect_identical( - dev_multinom(1, 10, 0.2), - dev_pois(1, 10 * 0.2) + dev_multinom(c(1, 4), c(10, 10), c(0.2, 0.8), group = c(1, 1)), + dev_pois(c(1, 4), c(10, 10) * c(0.2, 0.8)) ) # row-level deviance identity: summing the per-cell deviances within a # trial recovers the classic multinomial saturated-model deviance x <- c(1, 3, 6) size <- 10 prob <- c(0.2, 0.3, 0.5) + group <- c(1, 1, 1) expect_equal( - sum(dev_multinom(x, size, prob)), + sum(dev_multinom(x, size, prob, group)), 2 * sum(x * log(x / (size * prob))) ) + # group is validated: prob must sum to 1 per group, and every group must + # have at least 2 rows and match the modal row count + expect_error( + dev_multinom(4, 10, 1, group = 1), + "must contain at least 2 rows" + ) + expect_error( + dev_multinom(c(1, 2), c(10, 10), c(0.2, 0.2), group = c(1, 1)), + "`prob` must sum to 1" + ) }) test_that("dev_multinom matches glmnet deviance", { @@ -515,6 +526,7 @@ test_that("dev_multinom matches glmnet deviance", { })) colnames(y_wide) <- paste0("cat", seq_len(K)) + group <- rep(seq_len(n), each = K) for (lambda in c(0, 0.01, 0.05)) { fit <- glmnet::glmnet( x, @@ -529,7 +541,8 @@ test_that("dev_multinom matches glmnet deviance", { sum(dev_multinom( as.vector(t(y_wide)), size = 1, - prob = as.vector(t(prob_hat)) + prob = as.vector(t(prob_hat)), + group = group )), as.numeric(dev_glmnet) ) @@ -541,9 +554,10 @@ test_that("dev_multinom res", { x <- c(1, 3, 6) size <- 10 prob <- c(0.2, 0.3, 0.5) - res <- dev_multinom(x, size, prob, res = TRUE) + group <- c(1, 1, 1) + res <- dev_multinom(x, size, prob, group, res = TRUE) expect_equal(sign(res), sign(x - size * prob)) - expect_equal(sum(res^2), sum(dev_multinom(x, size, prob))) + expect_equal(sum(res^2), sum(dev_multinom(x, size, prob, group))) }) test_that("dev_neg_binom", { diff --git a/tests/testthat/test-res.R b/tests/testthat/test-res.R index 3d615875..4dffad57 100644 --- a/tests/testthat/test-res.R +++ b/tests/testthat/test-res.R @@ -515,7 +515,7 @@ test_that("res_multinom", { expect_equal( res_multinom(x, size, prob, group, type = "dev"), - dev_multinom(x, size, prob, res = TRUE) + dev_multinom(x, size, prob, group, res = TRUE) ) expect_equal( res_multinom(x, size, prob, group, type = "raw"), @@ -579,6 +579,16 @@ test_that("res_multinom", { ), "must not have any missing values" ) + # group is validated even when simulate = FALSE, so a bad prob sum can't + # silently fall through to a wrong (rather than an errored) residual + expect_error( + res_multinom(4, 10, 1, 1, simulate = FALSE), + "must contain at least 2 rows" + ) + expect_error( + res_multinom(c(4, 6), c(10, 10), c(0.4, 0.4), c(1, 1), simulate = FALSE), + "`prob` must sum to 1" + ) # sum of squared deviance residuals recovers the row-level deviance expect_equal( From 7e1663d51402ede50e9a7b44ef1708dc3fff856e Mon Sep 17 00:00:00 2001 From: Nicole Hill Date: Mon, 10 Aug 2026 09:58:01 -0700 Subject: [PATCH 7/9] Reword singleton-group error to hint at vectorized evaluation The most common cause is size/prob/group being passed one row at a time rather than as full vectors, so point at that instead of "group lost rows". --- R/internal.R | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/R/internal.R b/R/internal.R index 94788643..88caf22d 100644 --- a/R/internal.R +++ b/R/internal.R @@ -15,7 +15,7 @@ chk_multinom_group <- function(size, prob, group) { for (idx in groups) { if (length(idx) < 2L) { stop( - "Each `group` must contain at least 2 rows (a multinomial trial needs at least 2 categories); found a group with only 1 row. This usually means `group` lost the rows that should share a trial before reaching here.", + "Each `group` must contain at least 2 rows (a multinomial trial needs at least 2 categories); found a group with only 1 row. This usually means `group`/`size`/`prob` were passed one row at a time instead of as vectors.", call. = FALSE ) } From b5ad25d4227c64f45b5c070018c49e11fb662716 Mon Sep 17 00:00:00 2001 From: Nicole Hill Date: Mon, 10 Aug 2026 10:31:14 -0700 Subject: [PATCH 8/9] Fix length-mismatch and partial-NA validation gaps in multinomial functions dev_multinom()/log_lik_multinom()/res_multinom() now validate size/prob/ group lengths via chk_compatible_lengths() before recycling with rep_len(), matching ran_multinom() and skewnorm()/skewlnorm() -- previously a non-recyclable length mismatch (e.g. prob shorter than x) was silently mis-recycled by rep_len() instead of erroring, producing a wrong or misleading result. chk_multinom_group()'s prob-sum check no longer skips validation entirely whenever any prob in a group is NA; it now checks the known values don't already exceed 1, mirroring the size check's existing NA-tolerant logic. Also share one group-index split (multinom_split()) across chk_multinom_group()/multinom_row_na()/ran_multinom() instead of recomputing it up to three times per call, and add direct tests for the two internal helpers. --- R/dev.R | 1 + R/internal.R | 29 +++++++++++++---- R/log-lik.R | 9 +++--- R/ran.R | 12 +++---- R/res.R | 1 + tests/testthat/test-dev.R | 10 ++++++ tests/testthat/test-internal.R | 59 ++++++++++++++++++++++++++++++++++ tests/testthat/test-log-lik.R | 11 +++++++ tests/testthat/test-ran.R | 6 ++++ tests/testthat/test-res.R | 15 +++++++++ 10 files changed, 135 insertions(+), 18 deletions(-) create mode 100644 tests/testthat/test-internal.R diff --git a/R/dev.R b/R/dev.R index 38b2380e..d9fda82f 100644 --- a/R/dev.R +++ b/R/dev.R @@ -292,6 +292,7 @@ dev_lnorm <- function(x, meanlog = 0, sdlog = 1, res = FALSE) { #' @examples #' dev_multinom(c(1, 3, 6), size = 10, prob = c(0.2, 0.3, 0.5), group = c(1, 1, 1)) dev_multinom <- function(x, size = 1, prob, group, res = FALSE) { + chk_compatible_lengths(x, size, prob, group) n <- length(x) size <- rep_len(size, n) prob <- rep_len(prob, n) diff --git a/R/internal.R b/R/internal.R index 88caf22d..2c000a3f 100644 --- a/R/internal.R +++ b/R/internal.R @@ -2,6 +2,13 @@ dev_res <- function(x, mu, dev) { sign(x - mu) * sqrt(dev) } +# Row indices for each `group`, shared across the multinom_* helpers within +# a single call (chk_multinom_group(), multinom_row_na(), the sampling loop +# in ran_multinom()) so `group` isn't re-split by every one of them. +multinom_split <- function(group) { + split(seq_along(group), group) +} + # Checks every group shares one `size` and `prob` values summing to 1 # (required by rmultinom() and the deviance/log-lik identities), has >= 2 # rows (a trial needs >= 2 categories -- singletons usually mean `group` @@ -10,8 +17,7 @@ dev_res <- function(x, mu, dev) { # lost). Only non-NA values are compared, so lone NAs don't error here -- # see multinom_row_na(). Callers must chk_not_any_na(group) first; `group` # itself can't be NA-tolerant since it's what identifies the trial. -chk_multinom_group <- function(size, prob, group) { - groups <- split(seq_along(group), group) +chk_multinom_group <- function(size, prob, group, groups = multinom_split(group)) { for (idx in groups) { if (length(idx) < 2L) { stop( @@ -26,8 +32,17 @@ chk_multinom_group <- function(size, prob, group) { call. = FALSE ) } - prob_sum <- sum(prob[idx]) - if (!is.na(prob_sum) && abs(prob_sum - 1) > 1e-6) { + known_prob <- prob[idx][!is.na(prob[idx])] + known_prob_sum <- sum(known_prob) + # a group with a missing prob can only be validated one-sided: the known + # values must not already exceed 1, since a full sum-to-1 check would be + # (wrongly) skipped whenever any prob in the group is NA + prob_bad <- if (length(known_prob) < length(idx)) { + known_prob_sum > 1 + 1e-6 + } else { + abs(known_prob_sum - 1) > 1e-6 + } + if (prob_bad) { stop( "`prob` must sum to 1 for every `group` (multinomial trial).", call. = FALSE @@ -37,6 +52,8 @@ chk_multinom_group <- function(size, prob, group) { if (length(groups) > 1L) { group_sizes <- lengths(groups) size_counts <- table(group_sizes) + # ties are broken in favour of the smallest row count (table()'s names + # are sorted ascending, and which.max() takes the first maximum) mode_size <- as.integer(names(size_counts)[which.max(size_counts)]) bad <- group_sizes != mode_size if (any(bad)) { @@ -55,10 +72,10 @@ chk_multinom_group <- function(size, prob, group) { # Flags every row whose trial has an NA `size`/`prob` anywhere in the group, # since a trial's categories are scored/drawn jointly, not independently. -multinom_row_na <- function(size, prob, group) { +multinom_row_na <- function(size, prob, group, groups = multinom_split(group)) { bad <- is.na(size) | is.na(prob) result <- rep(FALSE, length(group)) - for (idx in split(seq_along(group), group)) { + for (idx in groups) { if (any(bad[idx])) { result[idx] <- TRUE } diff --git a/R/log-lik.R b/R/log-lik.R index cf88645e..044a61a0 100644 --- a/R/log-lik.R +++ b/R/log-lik.R @@ -457,22 +457,21 @@ log_lik_lnorm <- function(x, meanlog = 0, sdlog = 1, tlower = 0, tupper = Inf) { #' @examples #' log_lik_multinom(c(1, 3, 6), size = 10, prob = c(0.2, 0.3, 0.5), group = c(1, 1, 1)) log_lik_multinom <- function(x, size = 1, prob, group) { + chk_compatible_lengths(x, size, prob, group) n <- length(x) - if (!n || !length(size) || !length(prob) || !length(group)) { - return(numeric(0)) - } size <- rep_len(size, n) prob <- rep_len(prob, n) group <- rep_len(group, n) chk_not_any_na(group) - chk_multinom_group(size, prob, group) + groups <- multinom_split(group) + chk_multinom_group(size, prob, group, groups) mu <- size * prob log_lik <- log_lik_pois(x, mu) group_size <- table(group) k <- as.numeric(group_size[as.character(group)]) const <- log_lik_pois(size, size) log_lik <- log_lik - const / k - log_lik[multinom_row_na(size, prob, group)] <- NA_real_ + log_lik[multinom_row_na(size, prob, group, groups)] <- NA_real_ log_lik } diff --git a/R/ran.R b/R/ran.R index e6bb6a8f..0d5863ce 100644 --- a/R/ran.R +++ b/R/ran.R @@ -138,19 +138,17 @@ ran_lnorm <- function(n = 1, meanlog = 0, sdlog = 1) { #' @examples #' ran_multinom(size = 10, prob = c(0.2, 0.3, 0.5), group = c(1, 1, 1)) ran_multinom <- function(size = 1, prob, group) { + chk_compatible_lengths(size, prob, group) n <- length(prob) - if (!n) { - return(integer(0)) - } - chk_compatible_lengths(rep(1, n), size, group) size <- rep_len(size, n) prob <- rep_len(prob, n) group <- rep_len(group, n) chk_not_any_na(group) - chk_multinom_group(size, prob, group) - row_na <- multinom_row_na(size, prob, group) + groups <- multinom_split(group) + chk_multinom_group(size, prob, group, groups) + row_na <- multinom_row_na(size, prob, group, groups) x <- rep(NA_real_, n) - for (idx in split(seq_len(n), group)) { + for (idx in groups) { if (row_na[idx[1]]) { next } diff --git a/R/res.R b/R/res.R index a4f4f4a4..8b9e16fe 100644 --- a/R/res.R +++ b/R/res.R @@ -289,6 +289,7 @@ res_multinom <- function( simulate = FALSE ) { chk_string(type) + chk_compatible_lengths(x, size, prob, group) n <- length(x) size <- rep_len(size, n) prob <- rep_len(prob, n) diff --git a/tests/testthat/test-dev.R b/tests/testthat/test-dev.R index 784e5a81..b4fcb0cc 100644 --- a/tests/testthat/test-dev.R +++ b/tests/testthat/test-dev.R @@ -478,6 +478,16 @@ test_that("dev_multinom", { dev_multinom(numeric(0), numeric(0), numeric(0), group = numeric(0)), numeric(0) ) + # a mismatched, non-recyclable length errors clearly instead of being + # silently (and wrongly) recycled/NA-padded by rep_len() + expect_error( + dev_multinom(1:3, 10, numeric(0), group = c(1, 1, 1)), + "must be all zero length or the same length" + ) + expect_error( + dev_multinom(1:5, c(10, 10, 10), c(0.2, 0.3, 0.5), group = c(1, 1, 1, 2, 2)), + "must be all zero length or the same length" + ) expect_identical( dev_multinom(c(2, 3, 5), size = 10, prob = c(0.2, 0.3, 0.5), group = c(1, 1, 1)), c(0, 0, 0) diff --git a/tests/testthat/test-internal.R b/tests/testthat/test-internal.R new file mode 100644 index 00000000..1653a01c --- /dev/null +++ b/tests/testthat/test-internal.R @@ -0,0 +1,59 @@ +test_that("chk_multinom_group singleton group errors", { + expect_error( + chk_multinom_group(size = 10, prob = 1, group = 1), + "must contain at least 2 rows" + ) +}) + +test_that("chk_multinom_group size check ignores NA but catches known mismatches", { + expect_null(chk_multinom_group(c(10, 10, NA), c(0.2, 0.3, 0.5), c(1, 1, 1))) + expect_error( + chk_multinom_group(c(10, 20, NA), c(0.2, 0.3, 0.5), c(1, 1, 1)), + "`size` must be the same for every row belonging to the same `group`" + ) +}) + +test_that("chk_multinom_group prob check ignores NA but still catches known values that already break the sum", { + # a lone NA prob shouldn't skip validation of the other, already-known values + expect_error( + chk_multinom_group(c(10, 10, 10), c(0.9, 0.9, NA), c(1, 1, 1)), + "`prob` must sum to 1 for every `group`" + ) + # known values that don't yet exceed 1 are fine to leave for the NA to complete + expect_null(chk_multinom_group(c(10, 10, 10), c(0.4, 0.3, NA), c(1, 1, 1))) + # a fully-known group must still sum to exactly 1 + expect_null(chk_multinom_group(c(10, 10, 10), c(0.2, 0.3, 0.5), c(1, 1, 1))) + expect_error( + chk_multinom_group(c(10, 10, 10), c(0.2, 0.3, 0.4), c(1, 1, 1)), + "`prob` must sum to 1 for every `group`" + ) +}) + +test_that("chk_multinom_group modal row count check", { + expect_null( + chk_multinom_group( + c(10, 10, 10, 6, 6, 6), + c(0.2, 0.3, 0.5, 0.2, 0.3, 0.5), + c(1, 1, 1, 2, 2, 2) + ) + ) + expect_error( + chk_multinom_group( + c(10, 10, 10, 6, 6), + c(0.2, 0.3, 0.5, 0.5, 0.5), + c(1, 1, 1, 2, 2) + ), + "Every `group` should have the same number of rows" + ) +}) + +test_that("multinom_row_na flags a whole group when any size/prob in it is NA", { + expect_identical( + multinom_row_na(c(10, NA, 10, 10), c(0.2, 0.8, 0.4, 0.6), c(1, 1, 2, 2)), + c(TRUE, TRUE, FALSE, FALSE) + ) + expect_identical( + multinom_row_na(c(10, 10), c(0.2, 0.8), c(1, 1)), + c(FALSE, FALSE) + ) +}) diff --git a/tests/testthat/test-log-lik.R b/tests/testthat/test-log-lik.R index 90fb0038..37abfd1d 100644 --- a/tests/testthat/test-log-lik.R +++ b/tests/testthat/test-log-lik.R @@ -337,6 +337,17 @@ test_that("log_lik_multinom", { log_lik_multinom(numeric(0), numeric(0), numeric(0), numeric(0)), numeric(0) ) + # a mismatched, non-recyclable length errors clearly instead of being + # silently (and wrongly) recycled by rep_len(), or silently returning + # numeric(0) just because one argument happened to be empty + expect_error( + log_lik_multinom(1:3, 10, numeric(0), c(1, 1, 1)), + "must be all zero length or the same length" + ) + expect_error( + log_lik_multinom(1:5, c(10, 10, 10), c(0.2, 0.3, 0.5), c(1, 1, 1, 2, 2)), + "must be all zero length or the same length" + ) expect_error( log_lik_multinom(c(1, 3, 6), c(10, 10, 5), c(0.2, 0.3, 0.5), c(1, 1, 1)), "`size` must be the same for every row belonging to the same `group`" diff --git a/tests/testthat/test-ran.R b/tests/testthat/test-ran.R index 12afa41c..68b11318 100644 --- a/tests/testthat/test-ran.R +++ b/tests/testthat/test-ran.R @@ -112,6 +112,12 @@ test_that("ran_multinom", { ran_multinom(size = numeric(0), prob = numeric(0), group = numeric(0)), integer(0) ) + # a mismatched, non-recyclable length errors clearly instead of being + # silently (and wrongly) recycled by rep_len() + expect_error( + ran_multinom(size = c(10, 10, 10), prob = c(0.2, 0.3, 0.5), group = c(1, 1, 1, 2, 2)), + "must be all zero length or the same length" + ) expect_error( ran_multinom(size = c(10, 5), prob = c(0.5, 0.5), group = c(1, 1)), "`size` must be the same for every row belonging to the same `group`" diff --git a/tests/testthat/test-res.R b/tests/testthat/test-res.R index 4dffad57..0fbf9c5b 100644 --- a/tests/testthat/test-res.R +++ b/tests/testthat/test-res.R @@ -513,6 +513,21 @@ test_that("res_multinom", { prob <- c(0.2, 0.3, 0.5) group <- c(1, 1, 1) + expect_identical( + res_multinom(numeric(0), numeric(0), numeric(0), numeric(0)), + numeric(0) + ) + # a mismatched, non-recyclable length errors clearly instead of being + # silently (and wrongly) recycled/NA-padded by rep_len() + expect_error( + res_multinom(1:3, 10, numeric(0), c(1, 1, 1)), + "must be all zero length or the same length" + ) + expect_error( + res_multinom(1:5, c(10, 10, 10), c(0.2, 0.3, 0.5), c(1, 1, 1, 2, 2)), + "must be all zero length or the same length" + ) + expect_equal( res_multinom(x, size, prob, group, type = "dev"), dev_multinom(x, size, prob, group, res = TRUE) From 2569a14fd34aa703e43aaf41b9d3f80bab96709e Mon Sep 17 00:00:00 2001 From: Nicole Hill Date: Mon, 10 Aug 2026 10:57:31 -0700 Subject: [PATCH 9/9] Clarify internal test intent and pin down modal-size tie-break behavior Use expect_no_error() instead of expect_null() where the intent is "this call succeeds," not "this call returns NULL." Add a regression test for chk_multinom_group()'s documented tie-break rule (favours the smaller row count on an exact tie). --- tests/testthat/test-internal.R | 21 +++++++++++++++++---- 1 file changed, 17 insertions(+), 4 deletions(-) diff --git a/tests/testthat/test-internal.R b/tests/testthat/test-internal.R index 1653a01c..c6fa4034 100644 --- a/tests/testthat/test-internal.R +++ b/tests/testthat/test-internal.R @@ -6,7 +6,7 @@ test_that("chk_multinom_group singleton group errors", { }) test_that("chk_multinom_group size check ignores NA but catches known mismatches", { - expect_null(chk_multinom_group(c(10, 10, NA), c(0.2, 0.3, 0.5), c(1, 1, 1))) + expect_no_error(chk_multinom_group(c(10, 10, NA), c(0.2, 0.3, 0.5), c(1, 1, 1))) expect_error( chk_multinom_group(c(10, 20, NA), c(0.2, 0.3, 0.5), c(1, 1, 1)), "`size` must be the same for every row belonging to the same `group`" @@ -20,9 +20,9 @@ test_that("chk_multinom_group prob check ignores NA but still catches known valu "`prob` must sum to 1 for every `group`" ) # known values that don't yet exceed 1 are fine to leave for the NA to complete - expect_null(chk_multinom_group(c(10, 10, 10), c(0.4, 0.3, NA), c(1, 1, 1))) + expect_no_error(chk_multinom_group(c(10, 10, 10), c(0.4, 0.3, NA), c(1, 1, 1))) # a fully-known group must still sum to exactly 1 - expect_null(chk_multinom_group(c(10, 10, 10), c(0.2, 0.3, 0.5), c(1, 1, 1))) + expect_no_error(chk_multinom_group(c(10, 10, 10), c(0.2, 0.3, 0.5), c(1, 1, 1))) expect_error( chk_multinom_group(c(10, 10, 10), c(0.2, 0.3, 0.4), c(1, 1, 1)), "`prob` must sum to 1 for every `group`" @@ -30,7 +30,7 @@ test_that("chk_multinom_group prob check ignores NA but still catches known valu }) test_that("chk_multinom_group modal row count check", { - expect_null( + expect_no_error( chk_multinom_group( c(10, 10, 10, 6, 6, 6), c(0.2, 0.3, 0.5, 0.2, 0.3, 0.5), @@ -47,6 +47,19 @@ test_that("chk_multinom_group modal row count check", { ) }) +test_that("chk_multinom_group breaks an exact row-count tie toward the smaller size", { + # 2 groups of 2 rows, 2 groups of 3 rows -- documented to prefer the + # smaller count (table()'s ascending sort + which.max()'s first-match) + expect_error( + chk_multinom_group( + c(10, 10, 10, 10, 10, 10, 10, 10, 10, 10), + c(0.5, 0.5, 0.2, 0.3, 0.5, 0.5, 0.5, 0.2, 0.3, 0.5), + c(1, 1, 2, 2, 2, 3, 3, 4, 4, 4) + ), + "should have the same number of rows \\(2," + ) +}) + test_that("multinom_row_na flags a whole group when any size/prob in it is NA", { expect_identical( multinom_row_na(c(10, NA, 10, 10), c(0.2, 0.8, 0.4, 0.6), c(1, 1, 2, 2)),