Add multinomial distribution functions

DESCRIPTION

@@ -33,6 +33,7 @@ Suggests:
     covr,
     extraDistr,
     ggplot2,
+    glmnet,
     hms,
     knitr,
     memoise,

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)

R/dev.R

@@ -267,6 +267,41 @@ dev_lnorm <- function(x, meanlog = 0, sdlog = 1, res = FALSE) {
   dev_norm(log(x), mean = meanlog, sd = sdlog, res = res)
 }
 
+#' Multinomial Deviances
+#'
+#' 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 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), 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)
+  group <- rep_len(group, n)
+  chk_not_any_na(group)
+  chk_multinom_group(size, prob, group)
+  dev_pois(x, lambda = size * prob, res = res)
+}
+
 #' Negative Binomial Deviances
 #'
 #' @inheritParams params

R/internal.R

@@ -1,3 +1,84 @@
 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`
+# 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 = multinom_split(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`/`size`/`prob` were passed one row at a time instead of as vectors.",
+        call. = FALSE
+      )
+    }
+    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
+      )
+    }
+    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
+      )
+    }
+  }
+  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)) {
+      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
+      )
+    }
+  }
+}
+
+# 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, groups = multinom_split(group)) {
+  bad <- is.na(size) | is.na(prob)
+  result <- rep(FALSE, length(group))
+  for (idx in groups) {
+    if (any(bad[idx])) {
+      result[idx] <- TRUE
+    }
+  }
+  result
+}

R/log-lik.R

@@ -427,6 +427,54 @@ log_lik_lnorm <- function(x, meanlog = 0, sdlog = 1, tlower = 0, tupper = Inf) {
   log_lik
 }
 
+#' 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`. `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`.
+#' @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) {
+  chk_compatible_lengths(x, size, prob, group)
+  n <- length(x)
+  size <- rep_len(size, n)
+  prob <- rep_len(prob, n)
+  group <- rep_len(group, n)
+  chk_not_any_na(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, groups)] <- NA_real_
+  log_lik
+}
+
 #' Negative Binomial Log-Likelihood
 #'
 #' @inheritParams params

R/params.R

@@ -16,6 +16,11 @@
 #'   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 (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.

R/ran.R

@@ -113,6 +113,50 @@ ran_lnorm <- function(n = 1, meanlog = 0, sdlog = 1) {
   stats::rlnorm(n, meanlog = meanlog, sdlog = sdlog)
 }
 
+#' Multinomial Random Samples
+#'
+#' 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)`), 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`. `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
+#'
+#' @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)
+  size <- rep_len(size, n)
+  prob <- rep_len(prob, n)
+  group <- rep_len(group, n)
+  chk_not_any_na(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 groups) {
+    if (row_na[idx[1]]) {
+      next
+    }
+    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.

R/res.R

@@ -251,6 +251,65 @@ res_lnorm <- function(
   )
 }
 
+#' 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 (see
+#' [log_lik_multinom()] for details). `res_multinom()` returns one residual
+#' 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.
+#' @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)
+  chk_compatible_lengths(x, size, prob, group)
+  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)
+  }
+  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, group = group, res = TRUE),
+    chk_subset(x, c("data", "raw", "dev", "standardized"))
+  )
+}
+
 #' Negative Binomial Residuals
 #'
 #' @inheritParams params

_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`'

inst/WORDLIST

@@ -11,6 +11,7 @@ Numericise
 Numericize
 ORCID
 POSIXct
+Poissons
 Psychol
 Schaub
 Shachar