From 2b7915d3e08ff9f44fdf0dfc525a81d2288dd4f0 Mon Sep 17 00:00:00 2001
From: Josh Sumner <51797700+joshqsumner@users.noreply.github.com>
Date: Fri, 27 Feb 2026 09:28:00 -0600
Subject: [PATCH 01/25] lay out dirichlet helper with notes

---
 R/conj_dirichletHelpers.R | 66 +++++++++++++++++++++++++++++++++++++++
 1 file changed, 66 insertions(+)
 create mode 100644 R/conj_dirichletHelpers.R

diff --git a/R/conj_dirichletHelpers.R b/R/conj_dirichletHelpers.R
new file mode 100644
index 00000000..20718b1c
--- /dev/null
+++ b/R/conj_dirichletHelpers.R
@@ -0,0 +1,66 @@
+#' @description
+#' Internal function for calculating the dirichlet distribution underlying
+#' multinomially distributed single value traits (counts).
+#' @param s1 A named vector/list/data.frame of numerics drawn from a multinomial distribution.
+#' @examples
+#' out <- .conj_dirichlet_sv(
+#'   s1 = list("A" = 10, "B" = 10, "C" = 5),
+#'   cred.int.level = 0.95,
+#'   plot = TRUE
+#' )
+#' lapply(out, head)
+#' @keywords internal
+#' @noRd
+.conj_dirichlet_sv <- function(s1 = NULL, priors = NULL,
+                               support = NULL, cred.int.level = NULL,
+                               calculatingSupport = FALSE) {
+  #' `make default prior if none provided`
+  #' I think the default prior should be a flat dirichlet
+  #' lowest entropy would mean very low precision
+  #' interesting comments on jeffrey's prior here: https://arxiv.org/pdf/1504.02689
+  #` NOTE see examples 1.4 and 1.6 in linked paper for criticism of the jeffrey prior
+  #' in this context given large K and small N (which is probably an unlikely scenario
+  #' with the data that I'm designing this for, so we could just mention in the docs)
+  #' They actually come around to the explicit dirichlet example in ex 1.7 (although
+  #' on first reading I thought they were talking about dirichlet's the entire time?)
+  #' where `Alpha_i = 1/K for all i in K`.
+  #' I think this would be a fine prior for our settings, it's going to be very low
+  #' precision and flat. Besides, the use case we have in mind is:
+  #' Prior_0 is updated to Posterior_Before by Before-intervention data
+  #' Posterior_Before is updated to Posterior_After by After-Intervention data
+  #' Marginal Betas from Posterior_After and Posterior_Before are compared.
+  #'
+  #' NOTE as a sample input thing I don't think we are building for that scenario,
+  #' the way I imagine using that would be:
+  #' S1 = before_data, S2 = before_data + after_data
+  #' so that the S2 group is effectively updated for the before data, then again for
+  #' the after data.
+
+  #' `update dirichlet prior with sufficient statistics (counts)`
+
+  #' `Define support if it is missing`
+  #' support is the same as beta binomial
+
+  #' `Make Posterior draws`
+  #' Posterior draws used for ROPE comparisons mainly.
+  #' need to consider if it's better to draw from rdirichlet or rbeta
+
+  #' `calculate HDI`
+  #' note that this might change some since we'd have HDE/HDI for every
+  #' marginal beta
+
+  
+  #' `calculate HDE`
+  #' note that this might change some since we'd have HDE/HDI for every
+  #' marginal beta
+
+  
+  #' `Store summary`
+  #' note that this might change some since we'd have HDE/HDI for every
+  #' marginal beta
+
+  #' `Save s1 data for plotting`
+  #' note that this might change significantly from other distributions
+  #' because we want to test/plot the marginal betas to keep the dimensions
+  #' low enough to understand (and make the math tractable?)
+}

From 1e52f96aef27b01869bc0bd80fc65c15438a5659 Mon Sep 17 00:00:00 2001
From: Josh Sumner <51797700+joshqsumner@users.noreply.github.com>
Date: Fri, 27 Feb 2026 09:31:22 -0600
Subject: [PATCH 02/25] consistent function naming

---
 R/conj_dirichletHelpers.R | 5 ++++-
 1 file changed, 4 insertions(+), 1 deletion(-)

diff --git a/R/conj_dirichletHelpers.R b/R/conj_dirichletHelpers.R
index 20718b1c..7c98ded1 100644
--- a/R/conj_dirichletHelpers.R
+++ b/R/conj_dirichletHelpers.R
@@ -9,9 +9,12 @@
 #'   plot = TRUE
 #' )
 #' lapply(out, head)
+#' @details
+#' See Examples 1.4, 1.6, and 1.7 for thoughts on default dirichlet prior here
+#' https://arxiv.org/pdf/1504.02689
 #' @keywords internal
 #' @noRd
-.conj_dirichlet_sv <- function(s1 = NULL, priors = NULL,
+.conj_multinomial_sv <- function(s1 = NULL, priors = NULL,
                                support = NULL, cred.int.level = NULL,
                                calculatingSupport = FALSE) {
   #' `make default prior if none provided`

From c2fc67ce4824482fd28de5272d44a925054aa2e7 Mon Sep 17 00:00:00 2001
From: Josh Sumner <51797700+joshqsumner@users.noreply.github.com>
Date: Fri, 27 Feb 2026 09:34:12 -0600
Subject: [PATCH 03/25] naming functions consistently

---
 R/{conj_dirichletHelpers.R => conj_multinomialHelpers.R} | 4 +++-
 1 file changed, 3 insertions(+), 1 deletion(-)
 rename R/{conj_dirichletHelpers.R => conj_multinomialHelpers.R} (94%)

diff --git a/R/conj_dirichletHelpers.R b/R/conj_multinomialHelpers.R
similarity index 94%
rename from R/conj_dirichletHelpers.R
rename to R/conj_multinomialHelpers.R
index 7c98ded1..24f5a313 100644
--- a/R/conj_dirichletHelpers.R
+++ b/R/conj_multinomialHelpers.R
@@ -11,7 +11,9 @@
 #' lapply(out, head)
 #' @details
 #' See Examples 1.4, 1.6, and 1.7 for thoughts on default dirichlet prior here
-#' https://arxiv.org/pdf/1504.02689
+#' https://arxiv.org/pdf/1504.02689 , updating rule defined in
+#' The Compendium of Conjugate Priors (https://www.johndcook.com/CompendiumOfConjugatePriors.pdf)
+#' Section 5.1
 #' @keywords internal
 #' @noRd
 .conj_multinomial_sv <- function(s1 = NULL, priors = NULL,

From 592fec2e719f3726f1166bfc9d1c3ba8f45a55cb Mon Sep 17 00:00:00 2001
From: Josh Sumner <51797700+joshqsumner@users.noreply.github.com>
Date: Fri, 27 Feb 2026 09:50:55 -0600
Subject: [PATCH 04/25] notes on other functions that dir-multinom would need

---
 R/conj_multinomialHelpers.R | 24 ++++++++++++++++++++++++
 1 file changed, 24 insertions(+)

diff --git a/R/conj_multinomialHelpers.R b/R/conj_multinomialHelpers.R
index 24f5a313..75ca4700 100644
--- a/R/conj_multinomialHelpers.R
+++ b/R/conj_multinomialHelpers.R
@@ -50,6 +50,10 @@
   #' Posterior draws used for ROPE comparisons mainly.
   #' need to consider if it's better to draw from rdirichlet or rbeta
 
+  #' `Calculate density over support`
+  #' marginal beta densities, probably don't need all of them, just the ones
+  #' specified in the hypothesis
+
   #' `calculate HDI`
   #' note that this might change some since we'd have HDE/HDI for every
   #' marginal beta
@@ -69,3 +73,23 @@
   #' because we want to test/plot the marginal betas to keep the dimensions
   #' low enough to understand (and make the math tractable?)
 }
+
+#' Might define the plot helper functions here as well, I don't think I have anything else
+#' that is quite in this situation for how it's using marginals, that's why I have the
+#' bivariate options for some distributions but a shared bivariate plotting function.
+#' Overall example would be drawn from `conjugate_plot::.conj_general_plot` and might just
+#' be pre-processing.
+#' Still need to decide how those things are done, thinking that it would be specified in
+#' the hypothesis, so instead of `(s1) ">" (s2)` I'd have the input be `(s1)K_1 > (s2)K_3`?
+#' I think that if there are 2 samples then it would assume your hypothesis has to be between
+#' the samples (previous example) but if there is only 1 sample that you're comparing categories
+#' in that sample.
+#' 
+#' This will take some real review of how hypotheses are handled to implement.
+#' Currently the hypothesis is passed to `.post.prob.from.pdfs` and is just used for logic in picking
+#' which comparison to make between PDFs, so I could have the dirichlet helper parse the hypothesis
+#' string and only return PDFs of the marginal betas that are requested? I think then I could return
+#' two for each sample based on the hypothesis, then just always take the first from s1 and the
+#' second from s2? That would mean I'd need a helper higher upstream to parse that hypothesis since
+#' by the time the `.conj_distHelper` function sees data it only has 1 sample.
+#' Not 100% sure that works right now but that feels reasonable.

From 487ef364e21ae7ca39fc37fd48b7588e2289e41e Mon Sep 17 00:00:00 2001
From: Josh Sumner <51797700+joshqsumner@users.noreply.github.com>
Date: Mon, 2 Mar 2026 14:44:21 -0600
Subject: [PATCH 05/25] adding ellipses arguments to helpers, passing
 hypothesis to helpers for use with multinomial data

---
 R/conjugate.R                            |  2 +-
 R/conjugate_BernoulliHelpers.R           | 12 ++++++------
 R/conjugate_betaHelpers.R                | 20 ++++++++++----------
 R/conjugate_binomialHelpers.R            | 12 ++++++------
 R/conjugate_bivariate_lognormalHelpers.R |  2 +-
 R/conjugate_bivariate_uniformHelpers.R   |  4 ++--
 R/conjugate_exponentialHelpers.R         |  2 +-
 R/conjugate_gammaHelpers.R               |  3 ++-
 R/conjugate_gaussianHelpers.R            |  4 ++--
 R/conjugate_logNormal2Helpers.R          |  4 ++--
 R/conjugate_logNormalHelpers.R           |  7 ++-----
 R/conjugate_negBinHelpers.R              | 15 +++++----------
 R/conjugate_paretoHelpers.R              |  4 ++--
 R/conjugate_poissonHelpers.R             |  7 +------
 R/conjugate_tHelpers.R                   |  4 ++--
 R/conjugate_uniformHelpers.R             |  4 ++--
 R/conjugate_vonmises2Helpers.R           | 22 +++++++++++-----------
 R/conjugate_vonmisesHelpers.R            | 22 +++++++++++-----------
 18 files changed, 69 insertions(+), 81 deletions(-)

diff --git a/R/conjugate.R b/R/conjugate.R
index ac45eff4..7ab7a106 100644
--- a/R/conjugate.R
+++ b/R/conjugate.R
@@ -360,7 +360,7 @@ conjugate <- function(s1 = NULL, s2 = NULL,
       "bivariate_uniform", "bivariate_gaussian", "bivariate_lognormal"
     ))
     matched_fun <- get(paste0(".conj_", matched_arg, "_", vec_suffix))
-    res <- matched_fun(sample, prior, support, cred.int.level)
+    res <- matched_fun(sample, prior, support, cred.int.level, hypothesis)
     return(res)
   })
   #* `combine results into an object to return`
diff --git a/R/conjugate_BernoulliHelpers.R b/R/conjugate_BernoulliHelpers.R
index b381dac5..e4caacb6 100644
--- a/R/conjugate_BernoulliHelpers.R
+++ b/R/conjugate_BernoulliHelpers.R
@@ -13,7 +13,11 @@
 #' @noRd
 .conj_bernoulli_sv <- function(s1 = NULL, priors = NULL,
                                support = NULL, cred.int.level = NULL,
-                               calculatingSupport = FALSE) {
+                               calculatingSupport = FALSE, ...) {
+  #* `Define support if it is missing`
+  if (is.null(support) && calculatingSupport) {
+    return(c(0.0001, 0.9999))
+  }
   #* `make default prior if none provided`
   if (is.null(priors)) {
     priors <- list(a = 0.5, b = 0.5)
@@ -24,10 +28,6 @@
   #* `Update beta prior with sufficient statistics`
   a1_prime <- priors$a[1] + sum(s1)
   b1_prime <- priors$b[1] + sum(!s1)
-  #* `Define support if it is missing`
-  if (is.null(support) && calculatingSupport) {
-    return(c(0.0001, 0.9999))
-  }
   out <- list()
   #* `Make Posterior Draws`
   out$posteriorDraws <- rbeta(10000, a1_prime, b1_prime)
@@ -37,7 +37,7 @@
   out$pdf <- pdf1
   #* `calculate highest density interval`
   hdi1 <- qbeta(c((1 - cred.int.level) / 2, (1 - ((1 - cred.int.level) / 2))), a1_prime, b1_prime)
-  #* `calculate highest density estimate``
+  #* `calculate highest density estimate`
   hde1 <- .betaHDE(a1_prime, b1_prime)
   #* `Store summary`
   out$summary <- data.frame(HDE_1 = hde1, HDI_1_low = hdi1[1], HDI_1_high = hdi1[2])
diff --git a/R/conjugate_betaHelpers.R b/R/conjugate_betaHelpers.R
index d636dd4d..cb1cb8fb 100644
--- a/R/conjugate_betaHelpers.R
+++ b/R/conjugate_betaHelpers.R
@@ -21,15 +21,15 @@
 #' @noRd
 .conj_beta_mv <- function(s1 = NULL, priors = NULL,
                           support = NULL, cred.int.level = NULL,
-                          calculatingSupport = FALSE) {
-  #* `make default prior if none provided`
-  if (is.null(priors)) {
-    priors <- list(a = 0.5, b = 0.5)
-  }
+                          calculatingSupport = FALSE, ...) {
   #* `Define dense Support`
   if (is.null(support) && calculatingSupport) {
     return(c(0.0001, 0.9999))
   }
+  #* `make default prior if none provided`
+  if (is.null(priors)) {
+    priors <- list(a = 0.5, b = 0.5)
+  }
   out <- list()
   #* `Reorder columns if they are not in the numeric order`
   histColsBin <- as.numeric(sub("[a-zA-Z_.]+", "", colnames(s1)))
@@ -102,18 +102,18 @@
 #' @noRd
 .conj_beta_sv <- function(s1 = NULL, priors = NULL,
                           support = NULL, cred.int.level = NULL,
-                          calculatingSupport = FALSE) {
+                          calculatingSupport = FALSE, ...) {
   if (any(c(s1) > 1) || any(c(s1) < 0)) {
     stop("Beta Distribution is only defined on [0,1]")
   }
-  #* `make default prior if none provided`
-  if (is.null(priors)) {
-    priors <- list(a = 0.5, b = 0.5)
-  }
   #* `Define dense Support`
   if (is.null(support) && calculatingSupport) {
     return(c(0.0001, 0.9999))
   }
+  #* `make default prior if none provided`
+  if (is.null(priors)) {
+    priors <- list(a = 0.5, b = 0.5)
+  }
   out <- list()
 
   #* `get parameters for s1 using method of moments``
diff --git a/R/conjugate_binomialHelpers.R b/R/conjugate_binomialHelpers.R
index 899e3557..6ff9b237 100644
--- a/R/conjugate_binomialHelpers.R
+++ b/R/conjugate_binomialHelpers.R
@@ -13,7 +13,12 @@
 #' @noRd
 
 .conj_binomial_sv <- function(s1 = NULL, priors = NULL, support = NULL, cred.int.level = NULL,
-                              calculatingSupport = FALSE) {
+                              calculatingSupport = FALSE, ...) {
+  #* `Define dense Support`
+  #* `p parameter is beta distributed`
+  if (is.null(support) && calculatingSupport) {
+    return(c(0.0001, 0.9999))
+  }
   #* `check stopping conditions`
   s1 <- .conj_binomial_formatter(s1)
   #* `separate data into counts and trials`
@@ -29,11 +34,6 @@
   if (is.null(priors)) {
     priors <- list(a = 0.5, b = 0.5)
   }
-  #* `Define dense Support`
-  #* `p parameter is beta distributed`
-  if (is.null(support) && calculatingSupport) {
-    return(c(0.0001, 0.9999))
-  }
   out <- list()
   #* `Update priors with observed counts`
   a1_prime <- priors$a[1] + sum(s1_successes)
diff --git a/R/conjugate_bivariate_lognormalHelpers.R b/R/conjugate_bivariate_lognormalHelpers.R
index 4f2a87ac..57aff84a 100644
--- a/R/conjugate_bivariate_lognormalHelpers.R
+++ b/R/conjugate_bivariate_lognormalHelpers.R
@@ -12,7 +12,7 @@
 
 .conj_bivariate_lognormal_sv <- function(s1 = NULL, priors = NULL,
                                          support = NULL, cred.int.level = NULL,
-                                         calculatingSupport = FALSE) {
+                                         calculatingSupport = FALSE, ...) {
   out <- list()
   #* `make default prior if none provided`
   #* conjugate prior needs alpha, beta, mu, prec (or var or sd)
diff --git a/R/conjugate_bivariate_uniformHelpers.R b/R/conjugate_bivariate_uniformHelpers.R
index 372dea59..893f0773 100644
--- a/R/conjugate_bivariate_uniformHelpers.R
+++ b/R/conjugate_bivariate_uniformHelpers.R
@@ -15,7 +15,7 @@
 
 .conj_bivariate_uniform_sv <- function(s1 = NULL, priors = NULL,
                                        support = NULL, cred.int.level = NULL,
-                                       calculatingSupport = FALSE) {
+                                       calculatingSupport = FALSE, ...) {
   out <- list()
   #* `make default prior if none provided`
   #* conjugate prior needs r1, r2, and alpha
@@ -177,7 +177,7 @@
 
 .conj_bivariate_uniform_mv <- function(s1 = NULL, priors = NULL,
                                        support = NULL, cred.int.level = NULL,
-                                       calculatingSupport = FALSE) {
+                                       calculatingSupport = FALSE, ...) {
   out <- list()
   #* `make default prior if none provided`
   #* conjugate prior needs r1, r2, and alpha
diff --git a/R/conjugate_exponentialHelpers.R b/R/conjugate_exponentialHelpers.R
index a884c2ed..feb27020 100644
--- a/R/conjugate_exponentialHelpers.R
+++ b/R/conjugate_exponentialHelpers.R
@@ -12,7 +12,7 @@
 #' @noRd
 .conj_exponential_sv <- function(s1 = NULL, priors = NULL,
                                  support = NULL, cred.int.level = NULL,
-                                 calculatingSupport = FALSE) {
+                                 calculatingSupport = FALSE, ...) {
   out <- list()
   #* `make default prior if none provided`
   if (is.null(priors)) {
diff --git a/R/conjugate_gammaHelpers.R b/R/conjugate_gammaHelpers.R
index 662d561c..a7668c78 100644
--- a/R/conjugate_gammaHelpers.R
+++ b/R/conjugate_gammaHelpers.R
@@ -1,3 +1,4 @@
+
 #' @description
 #' Internal function for calculating the gamma distribution of the rate parameter in gamma distributed
 #' data represented by single value traits.
@@ -12,7 +13,7 @@
 #' @noRd
 .conj_gamma_sv <- function(s1 = NULL, priors = NULL,
                            support = NULL, cred.int.level = NULL,
-                           calculatingSupport = FALSE) {
+                           calculatingSupport = FALSE, ...) {
   out <- list()
   #* `make default prior if none provided`
   if (is.null(priors)) {
diff --git a/R/conjugate_gaussianHelpers.R b/R/conjugate_gaussianHelpers.R
index cc02a4d9..a2d62f3e 100644
--- a/R/conjugate_gaussianHelpers.R
+++ b/R/conjugate_gaussianHelpers.R
@@ -15,7 +15,7 @@
 
 .conj_gaussian_sv <- function(s1 = NULL, priors = NULL,
                               support = NULL, cred.int.level = NULL,
-                              calculatingSupport = FALSE) {
+                              calculatingSupport = FALSE, ...) {
   out <- list()
   #* `make default prior if none provided`
   priors <- .convert_gaussian_priors(priors)
@@ -76,7 +76,7 @@
 
 .conj_gaussian_mv <- function(s1 = NULL, priors = NULL,
                               support = NULL, cred.int.level = NULL,
-                              calculatingSupport = FALSE) {
+                              calculatingSupport = FALSE, ...) {
   out <- list()
   #* `make default prior if none provided`
   priors <- .convert_gaussian_priors(priors)
diff --git a/R/conjugate_logNormal2Helpers.R b/R/conjugate_logNormal2Helpers.R
index 25683b90..2ff42eb9 100644
--- a/R/conjugate_logNormal2Helpers.R
+++ b/R/conjugate_logNormal2Helpers.R
@@ -23,7 +23,7 @@
 
 .conj_lognormal2_mv <- function(s1 = NULL, priors = NULL,
                                 support = NULL, cred.int.level = NULL,
-                                calculatingSupport = FALSE) {
+                                calculatingSupport = FALSE, ...) {
   out <- list()
   #* `make default prior if none provided`
   if (is.null(priors)) {
@@ -110,7 +110,7 @@
 
 .conj_lognormal2_sv <- function(s1 = NULL, priors = NULL,
                                 support = NULL, cred.int.level = NULL,
-                                calculatingSupport = FALSE) {
+                                calculatingSupport = FALSE, ...) {
   out <- list()
   #* `make default prior if none provided`
   if (is.null(priors)) {
diff --git a/R/conjugate_logNormalHelpers.R b/R/conjugate_logNormalHelpers.R
index 56c503e5..944b6791 100644
--- a/R/conjugate_logNormalHelpers.R
+++ b/R/conjugate_logNormalHelpers.R
@@ -25,7 +25,7 @@
 
 .conj_lognormal_mv <- function(s1 = NULL, priors = NULL,
                                support = NULL, cred.int.level = NULL,
-                               calculatingSupport = FALSE) {
+                               calculatingSupport = FALSE, ...) {
   out <- list()
   #* `make default prior if none provided`
   priors <- .convert_gaussian_priors(priors)
@@ -36,9 +36,7 @@
   histColsBin <- as.numeric(sub("[a-zA-Z_.]+", "", colnames(s1)))
   bins_order <- sort(histColsBin, index.return = TRUE)$ix
   s1 <- s1[, bins_order]
-
   #* `Loop over reps, get moments for each histogram`
-
   rep_distributions <- lapply(seq_len(nrow(s1)), function(i) {
     X1 <- rep(histColsBin[bins_order], as.numeric(s1[i, ]))
     #* `Get mean of X1`
@@ -71,7 +69,6 @@
   obs_sd <- mean(unlist(lapply(rep_distributions, function(i) {
     return(1 / ((i$obs_prec) ^ 2))
   })))
-
   #* `Define support if it is missing`
   if (is.null(support) && calculatingSupport) {
     quantiles <- stats::qnorm(c(0.0001, 0.9999), mu_ls_prime, sigma_ls_prime)
@@ -129,7 +126,7 @@
 
 .conj_lognormal_sv <- function(s1 = NULL, priors = NULL,
                                support = NULL, cred.int.level = NULL,
-                               calculatingSupport = FALSE) {
+                               calculatingSupport = FALSE, ...) {
   out <- list()
   #* `make default prior if none provided`
   priors <- .convert_gaussian_priors(priors)
diff --git a/R/conjugate_negBinHelpers.R b/R/conjugate_negBinHelpers.R
index 11d1f521..adf07faf 100644
--- a/R/conjugate_negBinHelpers.R
+++ b/R/conjugate_negBinHelpers.R
@@ -31,7 +31,11 @@
 
 .conj_negbin_sv <- function(s1 = NULL, priors = NULL,
                             support = NULL, cred.int.level = NULL,
-                            calculatingSupport = FALSE) {
+                            calculatingSupport = FALSE, ...) {
+  #* `Define support if it is missing`
+  if (is.null(support) && calculatingSupport) {
+    return(c(0.0001, 0.9999))
+  }
   #* `Check samples`
   if (any(abs(s1 - round(s1)) > .Machine$double.eps^0.5) || any(s1 < 0)) {
     stop("Only positive whole numbers can be used in the Negative Binomial distribution")
@@ -44,27 +48,18 @@
       " you should add a prior including r parameter."
     ))
   }
-
   out <- list()
-
   #* `Use conjugate beta prior on probability`
   #* Note that this is very sensitive to the R value being appropriate
   a1_prime <- priors$a[1] + priors$r[1] * length(s1)
   b1_prime <- priors$b[1] + sum(s1)
-  #* `Define support if it is missing`
-  if (is.null(support) && calculatingSupport) {
-    return(c(0.0001, 0.9999))
-  }
   #* `calculate density over support``
   dens1 <- dbeta(support, a1_prime, b1_prime)
   pdf1 <- dens1 / sum(dens1)
-
   #* `calculate highest density interval`
   hdi1 <- qbeta(c((1 - cred.int.level) / 2, (1 - ((1 - cred.int.level) / 2))), a1_prime, b1_prime)
-
   #* `calculate highest density estimate``
   hde1 <- .betaHDE(a1_prime, b1_prime)
-
   #* `save summary and parameters`
   out$summary <- data.frame(HDE_1 = hde1, HDI_1_low = hdi1[1], HDI_1_high = hdi1[2])
   out$posterior$r <- priors$r[1]
diff --git a/R/conjugate_paretoHelpers.R b/R/conjugate_paretoHelpers.R
index 6fceb3c7..9f573462 100644
--- a/R/conjugate_paretoHelpers.R
+++ b/R/conjugate_paretoHelpers.R
@@ -19,7 +19,7 @@
 #' @noRd
 .conj_pareto_mv <- function(s1 = NULL, priors = NULL,
                             support = NULL, cred.int.level = NULL,
-                            calculatingSupport = FALSE) {
+                            calculatingSupport = FALSE, ...) {
   out <- list()
   #* `N observations`
   n_obs <- nrow(s1)
@@ -109,7 +109,7 @@
 #' @noRd
 .conj_pareto_sv <- function(s1 = NULL, priors = NULL,
                             support = NULL, cred.int.level = NULL,
-                            calculatingSupport = FALSE) {
+                            calculatingSupport = FALSE, ...) {
   out <- list()
   #* `make default prior if none provided`
   if (is.null(priors)) {
diff --git a/R/conjugate_poissonHelpers.R b/R/conjugate_poissonHelpers.R
index e48fe578..1011ad7f 100644
--- a/R/conjugate_poissonHelpers.R
+++ b/R/conjugate_poissonHelpers.R
@@ -23,7 +23,7 @@
 
 .conj_poisson_sv <- function(s1 = NULL, priors = NULL,
                              support = NULL, cred.int.level = NULL,
-                             calculatingSupport = FALSE) {
+                             calculatingSupport = FALSE, ...) {
   #* `Check samples`
   if (any(abs(s1 - round(s1)) > .Machine$double.eps^0.5) || any(s1 < 0)) {
     stop("Only positive integers can be used in the Poisson distribution")
@@ -32,9 +32,7 @@
   if (is.null(priors)) {
     priors <- list(a = 0.5, b = 0.5) # gamma prior on lambda
   }
-
   out <- list()
-
   #* `Use conjugate gamma prior on lambda`
   a1_prime <- priors$a[1] + sum(s1)
   b1_prime <- priors$b[1] + length(s1)
@@ -46,13 +44,10 @@
   #* `calculate density over support``
   dens1 <- dgamma(support, a1_prime, b1_prime)
   pdf1 <- dens1 / sum(dens1)
-
   #* `calculate highest density interval`
   hdi1 <- qgamma(c((1 - cred.int.level) / 2, (1 - ((1 - cred.int.level) / 2))), a1_prime, b1_prime)
-
   #* `calculate highest density estimate``
   hde1 <- .gammaHDE(shape = a1_prime, scale = 1 / b1_prime)
-
   #* `save summary and parameters`
   out$summary <- data.frame(HDE_1 = hde1, HDI_1_low = hdi1[1], HDI_1_high = hdi1[2])
   out$posterior$a <- a1_prime
diff --git a/R/conjugate_tHelpers.R b/R/conjugate_tHelpers.R
index ddd196bc..e4008781 100644
--- a/R/conjugate_tHelpers.R
+++ b/R/conjugate_tHelpers.R
@@ -14,7 +14,7 @@
 #' @noRd
 .conj_t_sv <- function(s1 = NULL, priors = NULL,
                        support = NULL, cred.int.level = NULL,
-                       calculatingSupport = FALSE) {
+                       calculatingSupport = FALSE, ...) {
   out <- list()
   #* `make default prior if none provided`
   priors <- .convert_gaussian_priors(priors)
@@ -93,7 +93,7 @@
 
 .conj_t_mv <- function(s1 = NULL, priors = NULL,
                        support = NULL, cred.int.level = NULL,
-                       calculatingSupport = FALSE) {
+                       calculatingSupport = FALSE, ...) {
   out <- list()
   #* `make default prior if none provided`
   priors <- .convert_gaussian_priors(priors)
diff --git a/R/conjugate_uniformHelpers.R b/R/conjugate_uniformHelpers.R
index cf71934e..d4c15a37 100644
--- a/R/conjugate_uniformHelpers.R
+++ b/R/conjugate_uniformHelpers.R
@@ -19,7 +19,7 @@
 #' @noRd
 .conj_uniform_mv <- function(s1 = NULL, priors = NULL,
                              support = NULL, cred.int.level = NULL,
-                             calculatingSupport = FALSE) {
+                             calculatingSupport = FALSE, ...) {
   out <- list()
   #* `make default prior if none provided`
   if (is.null(priors)) {
@@ -90,7 +90,7 @@
 #' @noRd
 .conj_uniform_sv <- function(s1 = NULL, priors = NULL,
                              support = NULL, cred.int.level = NULL,
-                             calculatingSupport = FALSE) {
+                             calculatingSupport = FALSE, ...) {
   out <- list()
   #* `make default prior if none provided`
   if (is.null(priors)) {
diff --git a/R/conjugate_vonmises2Helpers.R b/R/conjugate_vonmises2Helpers.R
index 28fd2889..410be41c 100644
--- a/R/conjugate_vonmises2Helpers.R
+++ b/R/conjugate_vonmises2Helpers.R
@@ -19,7 +19,7 @@
 
 .conj_vonmises2_sv <- function(s1 = NULL, priors = NULL,
                                support = NULL, cred.int.level = NULL,
-                               calculatingSupport = FALSE) {
+                               calculatingSupport = FALSE, ...) {
   #* `set support to NULL to avoid default length of 10000`
   support <- NULL
   #* `make default prior if none provided`
@@ -39,6 +39,15 @@
       return(default_prior[[nm]])
     }
   }), names(default_prior))
+  #* `Define dense Support`
+  if (is.null(support)) {
+    if (calculatingSupport) {
+      return(priors$boundary) #* this would be [-pi, pi] if using radians, but plotting will be on
+      #* the original scale so we can just return the boundary and use [-pi, pi] as support here
+    }
+    support_boundary <- seq(min(priors$boundary), max(priors$boundary), by = 0.0005)
+    support <- seq(-pi, pi, length.out = length(support_boundary))
+  }
   #* `rescale data to [-pi, pi] according to boundary`
   s1 <- .boundary.to.radians(x = s1, boundary = priors$boundary)
   #* `rescale prior on mu to [-pi, pi] according to boundary`
@@ -50,15 +59,6 @@
       "Does the boundary element in your prior include all your data?"
     ))
   }
-  #* `Define dense Support`
-  if (is.null(support)) {
-    if (calculatingSupport) {
-      return(priors$boundary) #* this would be [-pi, pi] if using radians, but plotting will be on
-      #* the original scale so we can just return the boundary and use [-pi, pi] as support here
-    }
-    support_boundary <- seq(min(priors$boundary), max(priors$boundary), by = 0.0005)
-    support <- seq(-pi, pi, length.out = length(support_boundary))
-  }
   out <- list()
   #* ***** `Updating Kappa`
   n1 <- length(s1)
@@ -146,7 +146,7 @@
 
 .conj_vonmises2_mv <- function(s1 = NULL, priors = NULL,
                                support = NULL, cred.int.level = NULL,
-                               calculatingSupport = FALSE) {
+                               calculatingSupport = FALSE, ...) {
   #* `set support to NULL to avoid default length of 10000`
   support <- NULL
   #* `Reorder columns if they are not in the numeric order`
diff --git a/R/conjugate_vonmisesHelpers.R b/R/conjugate_vonmisesHelpers.R
index 57787b1f..548ed191 100644
--- a/R/conjugate_vonmisesHelpers.R
+++ b/R/conjugate_vonmisesHelpers.R
@@ -20,7 +20,7 @@
 
 .conj_vonmises_mv <- function(s1 = NULL, priors = NULL,
                               support = NULL, cred.int.level = NULL,
-                              calculatingSupport = FALSE) {
+                              calculatingSupport = FALSE, ...) {
   #* `Turn off support for consistent rescaling between boundaries and to avoid default length of 10000`
   support <- NULL
   #* `Reorder columns if they are not in the numeric order`
@@ -46,6 +46,15 @@
       return(default_prior[[nm]])
     }
   }), names(default_prior))
+  #* `Define dense Support`
+  if (is.null(support)) {
+    if (calculatingSupport) {
+      return(priors$boundary) #* this would be [-pi, pi] if using radians, but plotting will be on
+      #* the original scale so we can just return the boundary and use [-pi, pi] as support here
+    }
+    support_boundary <- seq(min(priors$boundary), max(priors$boundary), by = 0.0005)
+    support <- seq(-pi, pi, length.out = length(support_boundary))
+  }
   #* `rescale data to [-pi, pi] according to boundary`
   X1 <- .boundary.to.radians(x = X1, boundary = priors$boundary)
   #* `rescale prior on mu to [-pi, pi] according to boundary`
@@ -57,15 +66,6 @@
       "Does the boundary element in your prior include all your data?"
     ))
   }
-  #* `Define dense Support`
-  if (is.null(support)) {
-    if (calculatingSupport) {
-      return(priors$boundary) #* this would be [-pi, pi] if using radians, but plotting will be on
-      #* the original scale so we can just return the boundary and use [-pi, pi] as support here
-    }
-    support_boundary <- seq(min(priors$boundary), max(priors$boundary), by = 0.0005)
-    support <- seq(-pi, pi, length.out = length(support_boundary))
-  }
   out <- list()
   #* `Get weighted mean of data and prior for half tangent adjustment`
   cm <- .circular.mean(c(X1, mu_radians), w = c(rep(nrow(s1) / length(X1), length(X1)), priors$n))
@@ -154,7 +154,7 @@
 
 .conj_vonmises_sv <- function(s1 = NULL, priors = NULL,
                               support = NULL, cred.int.level = NULL,
-                              calculatingSupport = FALSE) {
+                              calculatingSupport = FALSE, ...) {
   #* `to avoid default support length of 10000 which may not span boundary well`
   support <- NULL
   #* `make default prior if none provided`

From 5a990e96f7cf3f83b370d8b484753e2c7ed0dd45 Mon Sep 17 00:00:00 2001
From: Josh Sumner <51797700+joshqsumner@users.noreply.github.com>
Date: Mon, 2 Mar 2026 14:48:06 -0600
Subject: [PATCH 06/25] removing plot and support arguments fully

---
 R/conjugate.R | 11 ++---------
 1 file changed, 2 insertions(+), 9 deletions(-)

diff --git a/R/conjugate.R b/R/conjugate.R
index 7ab7a106..8cd6ced3 100644
--- a/R/conjugate.R
+++ b/R/conjugate.R
@@ -315,9 +315,9 @@ conjugate <- function(s1 = NULL, s2 = NULL,
                         "uniform", "pareto", "gamma", "bernoulli", "exponential", "bivariate_uniform",
                         "bivariate_gaussian", "bivariate_lognormal"
                       ),
-                      priors = NULL, plot = NULL, rope_range = NULL,
+                      priors = NULL, rope_range = NULL,
                       rope_ci = 0.89, cred.int.level = 0.89, hypothesis = "equal",
-                      bayes_factor = NULL, support = NULL) {
+                      bayes_factor = NULL) {
   #* `Handle formula option in s1`
   samples <- .formatSamples(s1, s2)
   s1 <- samples$s1
@@ -333,13 +333,6 @@ conjugate <- function(s1 = NULL, s2 = NULL,
   if (!is.null(s2)) {
     samplesList[[2]] <- s2
   }
-
-  if (!missing("support")) {
-    warning("support argument is deprecated.")
-  }
-  if (!missing("plot")) {
-    warning("plot argument is deprecated, use plot.conjugate instead.")
-  }
   support <- .getSupport(samplesList, method, priors) # calculate shared support
 
   sample_results <- lapply(seq_along(samplesList), function(i) {

From 094f9e8bdfd6110601089a0db428e2fc948e08ce Mon Sep 17 00:00:00 2001
From: Josh Sumner <51797700+joshqsumner@users.noreply.github.com>
Date: Wed, 15 Apr 2026 13:57:56 -0500
Subject: [PATCH 07/25] rename, add hypothesis helper

---
 ...lpers.R => conjugate_multinomialHelpers.R} | 123 ++++++++++++++++--
 1 file changed, 113 insertions(+), 10 deletions(-)
 rename R/{conj_multinomialHelpers.R => conjugate_multinomialHelpers.R} (54%)

diff --git a/R/conj_multinomialHelpers.R b/R/conjugate_multinomialHelpers.R
similarity index 54%
rename from R/conj_multinomialHelpers.R
rename to R/conjugate_multinomialHelpers.R
index 75ca4700..6095d721 100644
--- a/R/conj_multinomialHelpers.R
+++ b/R/conjugate_multinomialHelpers.R
@@ -1,7 +1,11 @@
 #' @description
 #' Internal function for calculating the dirichlet distribution underlying
 #' multinomially distributed single value traits (counts).
-#' @param s1 A named vector/list/data.frame of numerics drawn from a multinomial distribution.
+#'
+#' Note that this function returns posterior draws and densities from the marginal beta
+#' distributions, NOT from the dirichlet distribution.
+#'
+#' @param s1 A named vector/list of numerics drawn from a multinomial distribution. A data.frame will trigger the mv option, which will take column sums of the data then pass to the sv method.
 #' @examples
 #' out <- .conj_dirichlet_sv(
 #'   s1 = list("A" = 10, "B" = 10, "C" = 5),
@@ -9,6 +13,16 @@
 #'   plot = TRUE
 #' )
 #' lapply(out, head)
+#'
+#' s1 = list("A" = 10, "B" = 10, "C" = 5)
+#' priors = NULL
+#' support = seq(0.0001, 0.9999, length.out = 10000)
+#' calculatingSupport = FALSE
+#' cred.int.level = 0.95
+#' hypothesis = "A > B"
+#' sample_number = 1
+#'
+#' 
 #' @details
 #' See Examples 1.4, 1.6, and 1.7 for thoughts on default dirichlet prior here
 #' https://arxiv.org/pdf/1504.02689 , updating rule defined in
@@ -16,10 +30,24 @@
 #' Section 5.1
 #' @keywords internal
 #' @noRd
+
 .conj_multinomial_sv <- function(s1 = NULL, priors = NULL,
                                support = NULL, cred.int.level = NULL,
-                               calculatingSupport = FALSE) {
+                               calculatingSupport = FALSE, ...) {
+  nms <- names(s1)
+  s1 <- unlist(s1)
+
+  #' `Define support if it is missing`
+  #' support is the same as beta binomial
+  if (is.null(support) && calculatingSupport) {
+    return(c(0.0001, 0.9999))
+  }
+  #' `Check stopping conditions (similar to binomial)`
+  
   #' `make default prior if none provided`
+  if (is.null(priors)) {
+    priors <- list("alpha" = rep(1 / length(s1), length(s1)))
+  }
   #' I think the default prior should be a flat dirichlet
   #' lowest entropy would mean very low precision
   #' interesting comments on jeffrey's prior here: https://arxiv.org/pdf/1504.02689
@@ -42,36 +70,80 @@
   #' the after data.
 
   #' `update dirichlet prior with sufficient statistics (counts)`
-
-  #' `Define support if it is missing`
-  #' support is the same as beta binomial
+  alpha_prime <- priors$alpha + s1
+  
+  #' `Parse hypothesis to get marginal distributions`
+  #' We don't need to parse the hypothesis here since we'll return the entire set
+  #' of marginal posterior draws and densities.
+  
 
   #' `Make Posterior draws`
   #' Posterior draws used for ROPE comparisons mainly.
   #' need to consider if it's better to draw from rdirichlet or rbeta
+  #' I think I want betas
+  marginal_post <- do.call(cbind, lapply(seq_along(alpha_prime), function(i) {
+    return(rbeta(10000, alpha_prime[i], sum(alpha_prime[-i])))
+  }))
 
   #' `Calculate density over support`
   #' marginal beta densities, probably don't need all of them, just the ones
   #' specified in the hypothesis
+  marginal_dens <- do.call(cbind, lapply(seq_along(alpha_prime), function(i) {
+    d <- dbeta(support, alpha_prime[i], sum(alpha_prime[-i]))
+    return(d / sum(d))
+  }))
 
   #' `calculate HDI`
-  #' note that this might change some since we'd have HDE/HDI for every
+  #' note that this might change some since we could have HDE/HDI for every
   #' marginal beta
+  hdis <- do.call(cbind, lapply(seq_along(alpha_prime), function(i) {
+    hdi <- qbeta(c((1 - cred.int.level) / 2, (1 - ((1 - cred.int.level) / 2))),
+      alpha_prime[i], sum(alpha_prime[-i]))
+    return(hdi)
+  }))
 
-  
   #' `calculate HDE`
-  #' note that this might change some since we'd have HDE/HDI for every
+  #' note that this might change some since we could have HDE/HDI for every
   #' marginal beta
+  hdes <- do.call(cbind, lapply(seq_along(alpha_prime), function(i) {
+    hde <- .betaHDE(alpha_prime[i], sum(alpha_prime[-i]))
+    return(hde)
+  }))
+  colnames(hdes) <- colnames(hdis) <- colnames(marginal_dens) <- colnames(marginal_post) <- nms
 
-  
   #' `Store summary`
-  #' note that this might change some since we'd have HDE/HDI for every
+  #' note that this might change some since we could have HDE/HDI for every
   #' marginal beta
+  out <- list()
+  out$summary <- data.frame(
+    HDE_1 = as.numeric(hdes), HDI_1_low = hdis[1, ], HDI_1_high = hdis[2, ]
+  )
+  rownames(out$summary) <- nms
+  out$posterior$alpha <- alpha_prime
+  out$prior <- priors
+  out$posteriorDraws <- marginal_post
+  out$pdf <- marginal_dens
 
   #' `Save s1 data for plotting`
   #' note that this might change significantly from other distributions
   #' because we want to test/plot the marginal betas to keep the dimensions
   #' low enough to understand (and make the math tractable?)
+  out$plot_list <- list(
+    "range" = range(support),
+    "ddist_fun" = "stats::dbeta",
+    "priors" = list("shape1" = priors$alpha,
+      "shape2" = unlist(lapply(seq_along(priors$alpha), function(i) {
+        return(sum(priors$alpha[-i]))
+      }))
+    ),
+    "parameters" = list("shape1" = alpha_prime,
+      "shape2" = unlist(lapply(seq_along(alpha_prime), function(i) {
+        return(sum(alpha_prime[-i]))
+      }))
+    ),
+    "given" = NULL
+  )
+  return(out)
 }
 
 #' Might define the plot helper functions here as well, I don't think I have anything else
@@ -93,3 +165,34 @@
 #' second from s2? That would mean I'd need a helper higher upstream to parse that hypothesis since
 #' by the time the `.conj_distHelper` function sees data it only has 1 sample.
 #' Not 100% sure that works right now but that feels reasonable.
+
+#' @keywords internal
+#' @noRd
+.multinomial.pdf.handling <- function(sample_results, hypothesis) {
+  g1 <- trimws(gsub("(.*?)([>|<|=|!]{1,2})(.*)", "\\1", hypothesis))
+  g2 <- trimws(gsub("(.*?)([>|<|=|!]{1,2})(.*)", "\\3", hypothesis))
+  hyp <- trimws(gsub("(.*?)([>|<|=|!]{1,2})(.*)", "\\2", hypothesis))
+  hyp <- switch({hyp},
+    ">" = "greater", "<" = "lesser", "==" = "equal", "=" = "equal", "unequal"
+  )
+  if (length(sample_results) == 2) {
+    pdf.handling.output <- .post.prob.from.pdfs(
+      sample_results[[1]]$pdf[[g1]],
+      sample_results[[2]]$pdf[[g2]],
+      hyp
+    )
+  } else {
+    pdf.handling.output <- .post.prob.from.pdfs(
+      sample_results[[1]]$pdf[[g1]],
+      sample_results[[1]]$pdf[[g2]],
+      hyp
+    )
+  }
+  return(
+    list(
+      "pdf.handling.output" = pdf.handling.output,
+      "hyp" = hyp
+    )
+  )
+}
+

From 765ad38374ff327fcdda6348df04356c94d9f208 Mon Sep 17 00:00:00 2001
From: Josh Sumner <51797700+joshqsumner@users.noreply.github.com>
Date: Wed, 15 Apr 2026 13:58:11 -0500
Subject: [PATCH 08/25] add multinomial method

---
 R/conjugate.R | 25 +++++++++++++++++++------
 1 file changed, 19 insertions(+), 6 deletions(-)

diff --git a/R/conjugate.R b/R/conjugate.R
index 8cd6ced3..1a82cf2e 100644
--- a/R/conjugate.R
+++ b/R/conjugate.R
@@ -52,7 +52,11 @@
 #' in computing HDI for samples, defaults to 0.89.
 #' @param hypothesis Direction of a hypothesis if two samples are provided.
 #'  Options are "unequal", "equal", "greater", and "lesser",
-#'   read as "sample1 greater than sample2".
+#'   read as "sample1 greater than sample2". For the "multinomial" method the hypothesis
+#' should be specified as "Group1 >|<|==|!= Group2" and comparisons will be made using the marginal
+#' Beta distributions. If s2 is supplied then the hypothesis is read as
+#' "Group1 (from s1) >|<|==|!= Group2 (from s2)", if s2 is not supplied then both groups are taken
+#' from s1.
 #' @param bayes_factor Optional point or interval to evaluate bayes factors on. Note that this
 #' generally only makes sense to use if you have informative priors where the change in odds between
 #' prior and posterior is meaningful about the data. If this is non-NULL then columns of bayes factors
@@ -313,7 +317,7 @@ conjugate <- function(s1 = NULL, s2 = NULL,
                         "t", "gaussian", "beta", "binomial",
                         "lognormal", "lognormal2", "poisson", "negbin", "vonmises", "vonmises2",
                         "uniform", "pareto", "gamma", "bernoulli", "exponential", "bivariate_uniform",
-                        "bivariate_gaussian", "bivariate_lognormal"
+                        "bivariate_gaussian", "bivariate_lognormal", "multinomial"
                       ),
                       priors = NULL, rope_range = NULL,
                       rope_ci = 0.89, cred.int.level = 0.89, hypothesis = "equal",
@@ -350,10 +354,11 @@ conjugate <- function(s1 = NULL, s2 = NULL,
       "lognormal", "lognormal2", "poisson", "negbin",
       "vonmises", "vonmises2",
       "uniform", "pareto", "gamma", "bernoulli", "exponential",
-      "bivariate_uniform", "bivariate_gaussian", "bivariate_lognormal"
+      "bivariate_uniform", "bivariate_gaussian", "bivariate_lognormal",
+      "multinomial"
     ))
     matched_fun <- get(paste0(".conj_", matched_arg, "_", vec_suffix))
-    res <- matched_fun(sample, prior, support, cred.int.level, hypothesis)
+    res <- matched_fun(sample, prior, support, cred.int.level)
     return(res)
   })
   #* `combine results into an object to return`
@@ -370,7 +375,15 @@ conjugate <- function(s1 = NULL, s2 = NULL,
     }
     return(s)
   }))
-  if (!is.null(s2)) {
+  if (method[1] == "multinomial") {
+    # multinomial has special hypothesis handling, see conjugate_multinomialHelpers
+    mult_prob <- .multinomial.pdf.handling(sample_results, hypothesis)
+    out$summary <- cbind(
+      out$summary,
+      data.frame("hyp" = mult_prob$hyp,
+        "post.prob" = as.numeric(mult_prob$pdf.handling.output$$post.prob))
+    )
+  } else if (!is.null(s2)) {
     postProbRes <- .pdf.handling(sample_results[[1]]$pdf, sample_results[[2]]$pdf, hypothesis)
     out$summary <- cbind(
       out$summary,
@@ -620,7 +633,7 @@ conjugate <- function(s1 = NULL, s2 = NULL,
 #' @keywords internal
 #' @noRd
 
-.pdf.handling <- function(pdf1, pdf2, hypothesis) {
+.pdf.handling <- function(method, pdf1, pdf2, hypothesis) {
   if (is.list(pdf1) && is.list(pdf2)) {
     pdf.handling.output <- as.data.frame(do.call(rbind, lapply(seq_along(pdf1), function(i) {
       pdf <- .post.prob.from.pdfs(pdf1[[i]], pdf2[[i]], hypothesis)

From 2b95573a38c08b0b69919a09597c43d98f745f89 Mon Sep 17 00:00:00 2001
From: Josh Sumner <51797700+joshqsumner@users.noreply.github.com>
Date: Wed, 15 Apr 2026 13:58:19 -0500
Subject: [PATCH 09/25] allow additional arguments

---
 R/conjugate_bivariate_gaussianHelpers.R | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/R/conjugate_bivariate_gaussianHelpers.R b/R/conjugate_bivariate_gaussianHelpers.R
index 724241be..58a56cef 100644
--- a/R/conjugate_bivariate_gaussianHelpers.R
+++ b/R/conjugate_bivariate_gaussianHelpers.R
@@ -13,7 +13,7 @@
 
 .conj_bivariate_gaussian_sv <- function(s1 = NULL, priors = NULL,
                                         support = NULL, cred.int.level = NULL,
-                                        calculatingSupport = FALSE) {
+                                        calculatingSupport = FALSE, ...) {
   out <- list()
   #* `make default prior if none provided`
   #* conjugate prior needs alpha, beta, mu, prec (or var or sd)

From 64065015d1250edfc06134e7df7444bccc682256 Mon Sep 17 00:00:00 2001
From: Josh Sumner <51797700+joshqsumner@users.noreply.github.com>
Date: Wed, 15 Apr 2026 13:58:28 -0500
Subject: [PATCH 10/25] printing for multinomial

---
 R/conjugate_class.R | 3 ++-
 1 file changed, 2 insertions(+), 1 deletion(-)

diff --git a/R/conjugate_class.R b/R/conjugate_class.R
index 13faa9a8..9271e2b2 100644
--- a/R/conjugate_class.R
+++ b/R/conjugate_class.R
@@ -92,7 +92,8 @@ print.conjugatesummary <- function(x, ...) {
       exponential = list("Gamma", "Rate", "Exponential"),
       bivariate_uniform = list("Bivariate Pareto", "Boundaries", "Uniform"),
       bivariate_gaussian = list("Normal/Gamma", "Mu/Sd", "Normal"),
-      bivariate_lognormal = list("Normal/Gamma", "Mu/Sd", "Lognormal")
+      bivariate_lognormal = list("Normal/Gamma", "Mu/Sd", "Lognormal"),
+      multinomial = list("Dirichlet", "P (alpha)", "Multinomial")
     )
     method_statement <- paste0(
       method_list[[1]], # conjugate distribution

From 447c2c1dcde1667f5939490c5062694ddb96d39f Mon Sep 17 00:00:00 2001
From: Josh Sumner <51797700+joshqsumner@users.noreply.github.com>
Date: Wed, 15 Apr 2026 14:04:38 -0500
Subject: [PATCH 11/25] mv data option

---
 R/conjugate_multinomialHelpers.R | 17 ++++++++++++++++-
 1 file changed, 16 insertions(+), 1 deletion(-)

diff --git a/R/conjugate_multinomialHelpers.R b/R/conjugate_multinomialHelpers.R
index 6095d721..4365ce71 100644
--- a/R/conjugate_multinomialHelpers.R
+++ b/R/conjugate_multinomialHelpers.R
@@ -146,6 +146,22 @@
   return(out)
 }
 
+#' Multinomial "multi-value data", aka a dataframe of counts
+#' @param s1 A data frame of counts, columns should represent categories
+#' @keywords internal
+#' @noRd
+
+.conj_multinomial_mv <- function(s1 = NULL, priors = NULL,
+                               support = NULL, cred.int.level = NULL,
+                               calculatingSupport = FALSE, ...) {
+  nms <- colnames(s1)
+  s1 <- colSums(s1)
+  names(s1) <- nms
+  out <- .conj_multinomial_sv(s1, priors, support, cred.int.level, calculatingSupport)
+  return(out)
+}
+
+
 #' Might define the plot helper functions here as well, I don't think I have anything else
 #' that is quite in this situation for how it's using marginals, that's why I have the
 #' bivariate options for some distributions but a shared bivariate plotting function.
@@ -195,4 +211,3 @@
     )
   )
 }
-

From c29e38b6843ffd54a1e1c246e3f84756ea7f5064 Mon Sep 17 00:00:00 2001
From: Josh Sumner <51797700+joshqsumner@users.noreply.github.com>
Date: Wed, 15 Apr 2026 14:31:04 -0500
Subject: [PATCH 12/25] making conj object successfully

---
 R/conjugate.R                    | 6 +++---
 R/conjugate_multinomialHelpers.R | 8 ++++----
 2 files changed, 7 insertions(+), 7 deletions(-)

diff --git a/R/conjugate.R b/R/conjugate.R
index 1a82cf2e..62622967 100644
--- a/R/conjugate.R
+++ b/R/conjugate.R
@@ -375,15 +375,15 @@ conjugate <- function(s1 = NULL, s2 = NULL,
     }
     return(s)
   }))
-  if (method[1] == "multinomial") {
+  if (method[1] == "multinomial" && hypothesis != "equal") {
     # multinomial has special hypothesis handling, see conjugate_multinomialHelpers
     mult_prob <- .multinomial.pdf.handling(sample_results, hypothesis)
     out$summary <- cbind(
       out$summary,
       data.frame("hyp" = mult_prob$hyp,
-        "post.prob" = as.numeric(mult_prob$pdf.handling.output$$post.prob))
+        "post.prob" = as.numeric(mult_prob$pdf.handling.output$post.prob))
     )
-  } else if (!is.null(s2)) {
+  } else if (!is.null(s2) && method[1] != "multinomial") {
     postProbRes <- .pdf.handling(sample_results[[1]]$pdf, sample_results[[2]]$pdf, hypothesis)
     out$summary <- cbind(
       out$summary,
diff --git a/R/conjugate_multinomialHelpers.R b/R/conjugate_multinomialHelpers.R
index 4365ce71..a0dd16d3 100644
--- a/R/conjugate_multinomialHelpers.R
+++ b/R/conjugate_multinomialHelpers.R
@@ -193,14 +193,14 @@
   )
   if (length(sample_results) == 2) {
     pdf.handling.output <- .post.prob.from.pdfs(
-      sample_results[[1]]$pdf[[g1]],
-      sample_results[[2]]$pdf[[g2]],
+      sample_results[[1]]$pdf[, g1],
+      sample_results[[2]]$pdf[, g2],
       hyp
     )
   } else {
     pdf.handling.output <- .post.prob.from.pdfs(
-      sample_results[[1]]$pdf[[g1]],
-      sample_results[[1]]$pdf[[g2]],
+      sample_results[[1]]$pdf[, g1],
+      sample_results[[1]]$pdf[, g2],
       hyp
     )
   }

From bbe6017469729ba6a5fb72424145d9260c734a97 Mon Sep 17 00:00:00 2001
From: Josh Sumner <51797700+joshqsumner@users.noreply.github.com>
Date: Wed, 15 Apr 2026 14:54:23 -0500
Subject: [PATCH 13/25] pretty printing for alpha vector

---
 R/conjugate_class.R | 10 ++++++++--
 1 file changed, 8 insertions(+), 2 deletions(-)

diff --git a/R/conjugate_class.R b/R/conjugate_class.R
index 9271e2b2..53e6d7a9 100644
--- a/R/conjugate_class.R
+++ b/R/conjugate_class.R
@@ -112,7 +112,10 @@ print.conjugatesummary <- function(x, ...) {
       prior_statement <- paste0(
         "Sample ", i, " Prior ", method_list[[1]], "(",
         paste(
-          paste(names(prior), round(unlist(prior), 3), sep = " = "),
+          paste(names(prior),    
+            unlist(lapply(names(prior), function(nm) {
+              paste(round(prior[[nm]], 3), collapse = ", ")
+            })), sep = " = "),
           collapse = ", "
         ),
         ")\n"
@@ -120,7 +123,10 @@ print.conjugatesummary <- function(x, ...) {
       posterior_statement <- paste0(
         "\tPosterior ", method_list[[1]], "(",
         paste(
-          paste(names(post), round(unlist(post), 3), sep = " = "),
+          paste(names(post),    
+            unlist(lapply(names(post), function(nm) {
+              paste(round(post[[nm]], 3), collapse = ", ")
+            })), sep = " = "),
           collapse = ", "
         ),
         ")\n"

From c7a65f61bed82becd4a65d4658b74ae598beb75b Mon Sep 17 00:00:00 2001
From: Josh Sumner <51797700+joshqsumner@users.noreply.github.com>
Date: Wed, 15 Apr 2026 15:38:23 -0500
Subject: [PATCH 14/25] separating out helper for hypothesis parsing

---
 R/conjugate_multinomialHelpers.R | 13 ++++++++++++-
 1 file changed, 12 insertions(+), 1 deletion(-)

diff --git a/R/conjugate_multinomialHelpers.R b/R/conjugate_multinomialHelpers.R
index a0dd16d3..9a96fd7d 100644
--- a/R/conjugate_multinomialHelpers.R
+++ b/R/conjugate_multinomialHelpers.R
@@ -184,13 +184,24 @@
 
 #' @keywords internal
 #' @noRd
-.multinomial.pdf.handling <- function(sample_results, hypothesis) {
+.multinomial.parse.hypothesis <- function(hypothesis) {
   g1 <- trimws(gsub("(.*?)([>|<|=|!]{1,2})(.*)", "\\1", hypothesis))
   g2 <- trimws(gsub("(.*?)([>|<|=|!]{1,2})(.*)", "\\3", hypothesis))
   hyp <- trimws(gsub("(.*?)([>|<|=|!]{1,2})(.*)", "\\2", hypothesis))
   hyp <- switch({hyp},
     ">" = "greater", "<" = "lesser", "==" = "equal", "=" = "equal", "unequal"
   )
+  return(list(g1, g2, hyp))
+}
+
+
+#' @keywords internal
+#' @noRd
+.multinomial.pdf.handling <- function(sample_results, hypothesis) {
+  parsed <- .multinomial.parse.hypothesis(hypothesis)
+  g1 <- parsed[[1]]
+  g2 <- parsed[[2]]
+  hyp <- parsed[[3]]
   if (length(sample_results) == 2) {
     pdf.handling.output <- .post.prob.from.pdfs(
       sample_results[[1]]$pdf[, g1],

From 0f91776a6f6563f477350e623abdf2bf779b8518 Mon Sep 17 00:00:00 2001
From: Josh Sumner <51797700+joshqsumner@users.noreply.github.com>
Date: Thu, 16 Apr 2026 15:16:05 -0500
Subject: [PATCH 15/25] dirmultnom plotting and rope

---
 R/conjugate.R                    |  8 ++-
 R/conjugate_multinomialHelpers.R | 95 ++++++++++++++++++++++++++++++--
 2 files changed, 97 insertions(+), 6 deletions(-)

diff --git a/R/conjugate.R b/R/conjugate.R
index 62622967..354f8a23 100644
--- a/R/conjugate.R
+++ b/R/conjugate.R
@@ -392,7 +392,7 @@ conjugate <- function(s1 = NULL, s2 = NULL,
   }
   #* `parse output and do ROPE`
   if (!is.null(rope_range)) {
-    rope_res <- .conj_rope(sample_results, rope_range, rope_ci, method)
+    rope_res <- .conj_rope(sample_results, rope_range, rope_ci, method, hypothesis)
     out$summary <- cbind(out$summary, rope_res$summary)
     out$rope_df <- rope_res$rope_df
   }
@@ -543,7 +543,7 @@ conjugate <- function(s1 = NULL, s2 = NULL,
 #' @keywords internal
 #' @noRd
 .conj_rope <- function(sample_results, rope_range = c(-0.1, 0.1),
-                       rope_ci = 0.89, method) {
+                       rope_ci = 0.89, method, hypothesis) {
   #* `if bivariate then call the bivariate option`
   #* note this will return to .conj_rope but with a non-bivariate method
   if (any(grepl("bivariate", method))) {
@@ -553,6 +553,10 @@ conjugate <- function(s1 = NULL, s2 = NULL,
     )
     return(rope_res)
   }
+  #' `Format PDF from multinomial`
+  if (any(method == "multinomial")) {
+    sample_results <- .multinomial.rope.format(sample_results, hypothesis)
+  }
   #* `ROPE Comparison`
   rope_res <- list()
   if (!is.null(rope_range)) {
diff --git a/R/conjugate_multinomialHelpers.R b/R/conjugate_multinomialHelpers.R
index 9a96fd7d..b3e675c6 100644
--- a/R/conjugate_multinomialHelpers.R
+++ b/R/conjugate_multinomialHelpers.R
@@ -48,6 +48,7 @@
   if (is.null(priors)) {
     priors <- list("alpha" = rep(1 / length(s1), length(s1)))
   }
+  names(priors$alpha) <- nms
   #' I think the default prior should be a flat dirichlet
   #' lowest entropy would mean very low precision
   #' interesting comments on jeffrey's prior here: https://arxiv.org/pdf/1504.02689
@@ -132,14 +133,14 @@
     "range" = range(support),
     "ddist_fun" = "stats::dbeta",
     "priors" = list("shape1" = priors$alpha,
-      "shape2" = unlist(lapply(seq_along(priors$alpha), function(i) {
+      "shape2" = setNames(unlist(lapply(seq_along(priors$alpha), function(i) {
         return(sum(priors$alpha[-i]))
-      }))
+      })), nms)
     ),
     "parameters" = list("shape1" = alpha_prime,
-      "shape2" = unlist(lapply(seq_along(alpha_prime), function(i) {
+      "shape2" = setNames(unlist(lapply(seq_along(alpha_prime), function(i) {
         return(sum(alpha_prime[-i]))
-      }))
+      })), nms)
     ),
     "given" = NULL
   )
@@ -222,3 +223,89 @@
     )
   )
 }
+
+
+#' @keywords internal
+#' @noRd
+.multinomial.conj.plot.format <- function(res) {
+#' not immediately sure how this should work yet with the hypothesis being shuttled around
+  #' basically what I want to have happen is that this selects only the PDF/posteriors from the
+  #' groups that are in the hypothesis and plots those as betas. How to pass the hypothesis
+  #' information is tricky. Could label it earlier on? Could have an extra attribute in `res`
+  #' to store the formatted hypothesis?
+  hypl <- .multinomial.parse.hypothesis(res$call$hypothesis)
+  g1 <- hypl[[1]]
+  g2 <- hypl[[2]]
+  hyp <- hypl[[3]]
+  new_res <- res
+  if (length(res$plot_parameters) == 2) {
+    # keep only relevant marginal beta parameters
+    new_res$plot_parameters[[1]]$parameters$shape1 <- res$plot_parameters[[1]]$parameters$shape1[g1]
+    new_res$plot_parameters[[2]]$parameters$shape1 <- res$plot_parameters[[2]]$parameters$shape1[g2]
+    new_res$plot_parameters[[1]]$parameters$shape2 <- res$plot_parameters[[1]]$parameters$shape2[g1]
+    new_res$plot_parameters[[2]]$parameters$shape2 <- res$plot_parameters[[2]]$parameters$shape2[g2]
+    # keep only relevant marginal priors
+    new_res$plot_parameters[[1]]$priors$shape1 <- res$plot_parameters[[1]]$priors$shape1[g1]
+    new_res$plot_parameters[[2]]$priors$shape1 <- res$plot_parameters[[2]]$priors$shape1[g2]
+    new_res$plot_parameters[[1]]$priors$shape2 <- res$plot_parameters[[1]]$priors$shape2[g1]
+    new_res$plot_parameters[[2]]$priors$shape2 <- res$plot_parameters[[2]]$priors$shape2[g1]
+    # subset summary to draw HDE/HDI correctly
+    new_res$summary <- cbind(
+      res$summary[g1, grepl("_1", colnames(res$summary))],
+      res$summary[g2, grepl("_2", colnames(res$summary))],
+      res$summary[1, !grepl("[1|2]", colnames(res$summary))]
+    )
+  } else {
+    # keep only relevant marginal beta parameters
+    new_res$plot_parameters[[1]]$parameters$shape1 <- res$plot_parameters[[1]]$parameters$shape1[g1]
+    new_res$plot_parameters[[1]]$parameters$shape1 <- res$plot_parameters[[1]]$parameters$shape2[g2]
+    # keep only relevant marginal priors
+    new_res$plot_parameters[[1]]$priors$shape1 <- res$plot_parameters[[1]]$priors$shape1[g1]
+    new_res$plot_parameters[[1]]$priors$shape1 <- res$plot_parameters[[1]]$priors$shape2[g2]
+    # subset summary to draw HDE/HDI correctly
+    ls()
+    g1_cols <- res$summary[g1, grepl("_1", colnames(res$summary))]
+    g2_cols <- res$summary[g2, grepl("_1", colnames(res$summary))]
+    colnames(g2_cols) <- gsub("_1", "_2", colnames(g2_cols))
+    new_res$summary <- cbind(
+      g1_cols,
+      g2_cols,
+      res$summary[1, !grepl("[1|2]", colnames(res$summary))]
+    )
+  }
+  return(new_res)
+}
+
+
+#' @keywords internal
+#' @noRd
+
+.multinomial.rope.format <- function(sample_results, hypothesis) {
+  hypl <- .multinomial.parse.hypothesis(hypothesis)
+  g1 <- hypl[[1]]
+  g2 <- hypl[[2]]
+  post1 <- sample_results[[1]]$posteriorDraws[, g1]
+  #' for ROPE comparison everything expects 2L results, so force that.
+  if (length(sample_results) == 1) {
+    post2 <- sample_results[[1]]$posteriorDraws[, g2]
+    new_sample_results <- list(sample_results, sample_results)
+  } else {
+    post2 <- sample_results[[2]]$posteriorDraws[, g2]
+    new_sample_results <- sample_results
+  }
+  new_sample_results[[1]]$posteriorDraws <- post1
+  new_sample_results[[2]]$posteriorDraws <- post2
+  return(new_sample_results)
+}
+
+#' @param s_res results from conjugate function thus far, currently the plot_list
+#' (for the distribution function name and values) element is all that is used.
+#' Internally this object is called `sample_results` in conjugate and only has
+#' one sample at a time passed to this function. Here it just needs to be
+#' formatted to only have the relevant pieces.
+#' @keywords internal
+#' @noRd
+
+.multinomial.bayes.factor.format <- function(s_res) {
+
+}

From dd9c17c27dc32548c0c6b5b1ec2a031a73c85c95 Mon Sep 17 00:00:00 2001
From: Josh Sumner <51797700+joshqsumner@users.noreply.github.com>
Date: Thu, 16 Apr 2026 15:16:19 -0500
Subject: [PATCH 16/25] apply format to plot

---
 R/conjugate_plot.R | 3 +++
 1 file changed, 3 insertions(+)

diff --git a/R/conjugate_plot.R b/R/conjugate_plot.R
index 4aed592f..c6e9b364 100644
--- a/R/conjugate_plot.R
+++ b/R/conjugate_plot.R
@@ -54,6 +54,9 @@ plot.conjugate <- function(x, ...) {
   if (any(grepl("bivariate", method))) {
     p <- .conj_bivariate_plot(res, rope_df, rope_range, rope_ci, dirSymbol)
   } else {
+    if (any(grepl("multinomial", method))) {
+      res <- .multinomial.conj.plot.format(res)
+    }
     p <- .conj_general_plot(res, rope_df, rope_range, rope_ci, dirSymbol, support, bayes_factor)
   }
   return(p)

From ac43118a923c923129abda4de6d4b75528406cb4 Mon Sep 17 00:00:00 2001
From: Josh Sumner <51797700+joshqsumner@users.noreply.github.com>
Date: Thu, 23 Apr 2026 13:23:58 -0500
Subject: [PATCH 17/25] changing bayes factor helper arguments

---
 R/conjugate_bayes_factors.R | 13 ++++++++-----
 1 file changed, 8 insertions(+), 5 deletions(-)

diff --git a/R/conjugate_bayes_factors.R b/R/conjugate_bayes_factors.R
index 3375097f..d028ef26 100644
--- a/R/conjugate_bayes_factors.R
+++ b/R/conjugate_bayes_factors.R
@@ -4,10 +4,12 @@
 #' Function to calcualte Bayes Factors using single or multi value traits with
 #' several distributions in the conjugate function.
 #' @param bayes_factor bayes factor range/point hypothesis passed from conjugate
-#' @param s_res results from conjugate function thus far, currently the plot_list
-#' (for the distribution function name and values) element is all that is used.
-#' Internally this object is called `sample_results` in conjugate and only has
-#' one sample at a time passed to this function.
+#' @param sample_results results from conjugate function thus far,
+#' currently the `plot_list` (for the distribution function name and values)
+#' element is all that is used. Internally this object is called
+#' `sample_results` in conjugate and only has
+#' one sample at a time used in this function.
+#' @param i numeric, which element of sample_results to use.
 #' @examples
 #' sample_results <- list(
 #' # other things that we don't need to use for this function... ,
@@ -27,7 +29,8 @@
 #' @keywords internal
 #' @noRd
 
-.conj_bayes_factor <- function(bayes_factor, s_res) {
+.conj_bayes_factor <- function(bayes_factor, sample_results, i) {
+  s_res <- sample_results[[i]]
   if (length(bayes_factor) == 1) {
     # point hypothesis
     post_args <- append(bayes_factor, s_res$plot_list$parameters)

From 0ef784e8a71aa64a40cb7f4a308b215ed57f44ec Mon Sep 17 00:00:00 2001
From: Josh Sumner <51797700+joshqsumner@users.noreply.github.com>
Date: Thu, 23 Apr 2026 13:28:02 -0500
Subject: [PATCH 18/25] docs update

---
 R/conjugate.R                    |  43 ++++++++++--
 R/conjugate_multinomialHelpers.R | 117 ++++++++-----------------------
 man/conjugate.Rd                 |  53 +++++++++++---
 3 files changed, 106 insertions(+), 107 deletions(-)

diff --git a/R/conjugate.R b/R/conjugate.R
index 354f8a23..7dfc87ba 100644
--- a/R/conjugate.R
+++ b/R/conjugate.R
@@ -15,9 +15,9 @@
 #' @param s2 An optional second sample, or if s1 is a formula then this should be a dataframe.
 #' This sample is shown in blue if plotted.
 #' @param method The distribution/method to use.
-#' Currently "t", "gaussian", "beta", "binomial", "lognormal", "lognormal2", "poisson",
-#' "negbin" (negative binomial), "uniform", "pareto", "gamma", "bernoulli", "exponential",
-#' "vonmises", and "vonmises2" are supported.
+#' Currently "bernoulli", "beta", "binomial", "exponential", "gamma", "gaussian",
+#' "lognormal", "lognormal2", "multinomial", "negbin" (negative binomial), "pareto",
+#' "poisson", "t", "uniform", "vonmises", and "vonmises2" are supported.
 #' The count (binomial, poisson and negative binomial), bernoulli, exponential,
 #' and pareto distributions are only implemented for single value traits due to their updating
 #' and/or the nature of the input data.
@@ -56,7 +56,9 @@
 #' should be specified as "Group1 >|<|==|!= Group2" and comparisons will be made using the marginal
 #' Beta distributions. If s2 is supplied then the hypothesis is read as
 #' "Group1 (from s1) >|<|==|!= Group2 (from s2)", if s2 is not supplied then both groups are taken
-#' from s1.
+#' from s1. For the multinomial method groups should be specified, so the hypothesis is written as
+#' "group1 > group2", where "group1" and "group2" would be informed by s1 unless s2 is provided in
+#' which case "group1" would be from s1 and "group2" would be from s2.
 #' @param bayes_factor Optional point or interval to evaluate bayes factors on. Note that this
 #' generally only makes sense to use if you have informative priors where the change in odds between
 #' prior and posterior is meaningful about the data. If this is non-NULL then columns of bayes factors
@@ -131,6 +133,12 @@
 #'     kappa from the data and updates the kappa prior as a weighted average of the data and the prior.
 #'     The mu parameter is then updated per Von-Mises conjugacy.
 #'     }
+#'     \item{\strong{"multinomial": }
+#'     \code{list(alpha = list("alpha" = rep(1/N_groups, N_groups))}, where alpha is the concentration
+#'     vector of the conjugate dirichlet distribution. For the multinomial method hypotheses are
+#'     specified with group names, so instead of "equal" the hypothesis could be
+#'     "genotypeX == genotypeY".
+#'     }
 #'     \item{\strong{"bivariate_uniform": }
 #'     \code{list(location_l = 1, location_u = 2, scale = 1)}, where scale is the
 #'     shared scale parameter of the pareto distributed upper and lower boundaries and location l and u
@@ -260,6 +268,26 @@
 #'   cred.int.level = 0.89, hypothesis = "equal"
 #' )
 #'
+#' # dirichlet-multinomial sv example
+#'
+#' dm_sv_ex <- conjugate(
+#'   s1 = list("A" = 10, "B" = 10, "C" = 5),
+#'   s2 = list("A" = 4, "B" = 12, "C" = 9),
+#'   method = "multinomial",
+#'   hypothesis = "A > A",
+#'   rope_range = c(-0.1, 0.1)
+#' )
+#'
+#' # dirichlet-multinomial mv example
+#'
+#' dm_mv_ex <- conjugate(
+#'   s1 = data.frame("A" = c(5,5), "B" = c(5, 5), "C" = c(2,3)),
+#'   s2 = data.frame("A" = c(2,2), "B" = c(7, 5), "C" = c(8,1)),
+#'   method = "multinomial",
+#'   hypothesis = "A > A",
+#'   rope_range = c(-0.1, 0.1)
+#' )
+#'
 #' # von mises mv example
 #'
 #' mv_gauss <- mvSim(
@@ -366,7 +394,7 @@ conjugate <- function(s1 = NULL, s2 = NULL,
   out$summary <- do.call(cbind, lapply(seq_along(sample_results), function(i) {
     s <- sample_results[[i]]$summary
     if (!is.null(bayes_factor)) { #* `Calculate Bayes Factors`
-      bf <- .conj_bayes_factor(bayes_factor, sample_results[[i]])
+      bf <- .conj_bayes_factor(bayes_factor, sample_results, i)
       s$bf_1 <- bf
     }
     if (i == 2) {
@@ -375,7 +403,7 @@ conjugate <- function(s1 = NULL, s2 = NULL,
     }
     return(s)
   }))
-  if (method[1] == "multinomial" && hypothesis != "equal") {
+  if (method[1] == "multinomial" && hypothesis != "equal") { # if hypothesis is the default then skip
     # multinomial has special hypothesis handling, see conjugate_multinomialHelpers
     mult_prob <- .multinomial.pdf.handling(sample_results, hypothesis)
     out$summary <- cbind(
@@ -542,6 +570,7 @@ conjugate <- function(s1 = NULL, s2 = NULL,
 #' this should take outputs from conjHelpers and compare the $posteriorDraws.
 #' @keywords internal
 #' @noRd
+
 .conj_rope <- function(sample_results, rope_range = c(-0.1, 0.1),
                        rope_ci = 0.89, method, hypothesis) {
   #* `if bivariate then call the bivariate option`
@@ -553,7 +582,7 @@ conjugate <- function(s1 = NULL, s2 = NULL,
     )
     return(rope_res)
   }
-  #' `Format PDF from multinomial`
+  #* `Format PDF from multinomial`
   if (any(method == "multinomial")) {
     sample_results <- .multinomial.rope.format(sample_results, hypothesis)
   }
diff --git a/R/conjugate_multinomialHelpers.R b/R/conjugate_multinomialHelpers.R
index b3e675c6..dc6a5c09 100644
--- a/R/conjugate_multinomialHelpers.R
+++ b/R/conjugate_multinomialHelpers.R
@@ -37,84 +37,58 @@
   nms <- names(s1)
   s1 <- unlist(s1)
 
-  #' `Define support if it is missing`
-  #' support is the same as beta binomial
+  #* `Define support if it is missing`
+  #* support is the same as beta binomial
   if (is.null(support) && calculatingSupport) {
     return(c(0.0001, 0.9999))
   }
-  #' `Check stopping conditions (similar to binomial)`
+  #* `Check stopping conditions (similar to binomial)`
   
-  #' `make default prior if none provided`
+  #* `make default prior if none provided`
   if (is.null(priors)) {
     priors <- list("alpha" = rep(1 / length(s1), length(s1)))
   }
   names(priors$alpha) <- nms
-  #' I think the default prior should be a flat dirichlet
-  #' lowest entropy would mean very low precision
-  #' interesting comments on jeffrey's prior here: https://arxiv.org/pdf/1504.02689
-  #` NOTE see examples 1.4 and 1.6 in linked paper for criticism of the jeffrey prior
-  #' in this context given large K and small N (which is probably an unlikely scenario
-  #' with the data that I'm designing this for, so we could just mention in the docs)
-  #' They actually come around to the explicit dirichlet example in ex 1.7 (although
-  #' on first reading I thought they were talking about dirichlet's the entire time?)
-  #' where `Alpha_i = 1/K for all i in K`.
-  #' I think this would be a fine prior for our settings, it's going to be very low
-  #' precision and flat. Besides, the use case we have in mind is:
-  #' Prior_0 is updated to Posterior_Before by Before-intervention data
-  #' Posterior_Before is updated to Posterior_After by After-Intervention data
-  #' Marginal Betas from Posterior_After and Posterior_Before are compared.
-  #'
-  #' NOTE as a sample input thing I don't think we are building for that scenario,
-  #' the way I imagine using that would be:
-  #' S1 = before_data, S2 = before_data + after_data
-  #' so that the S2 group is effectively updated for the before data, then again for
-  #' the after data.
 
-  #' `update dirichlet prior with sufficient statistics (counts)`
+  #* `update dirichlet prior with sufficient statistics (counts)`
   alpha_prime <- priors$alpha + s1
-  
-  #' `Parse hypothesis to get marginal distributions`
-  #' We don't need to parse the hypothesis here since we'll return the entire set
-  #' of marginal posterior draws and densities.
-  
-
-  #' `Make Posterior draws`
-  #' Posterior draws used for ROPE comparisons mainly.
-  #' need to consider if it's better to draw from rdirichlet or rbeta
-  #' I think I want betas
+  #* `Make Posterior draws`
+  #* Posterior draws used for ROPE comparisons mainly.
+  #* need to consider if it's better to draw from rdirichlet or rbeta
+  #* I think I want betas
   marginal_post <- do.call(cbind, lapply(seq_along(alpha_prime), function(i) {
     return(rbeta(10000, alpha_prime[i], sum(alpha_prime[-i])))
   }))
 
-  #' `Calculate density over support`
-  #' marginal beta densities, probably don't need all of them, just the ones
-  #' specified in the hypothesis
+  #* `Calculate density over support`
+  #* marginal beta densities, probably don't need all of them, just the ones
+  #* specified in the hypothesis
   marginal_dens <- do.call(cbind, lapply(seq_along(alpha_prime), function(i) {
     d <- dbeta(support, alpha_prime[i], sum(alpha_prime[-i]))
     return(d / sum(d))
   }))
 
-  #' `calculate HDI`
-  #' note that this might change some since we could have HDE/HDI for every
-  #' marginal beta
+  #* `calculate HDI`
+  #* note that this might change some since we could have HDE/HDI for every
+  #* marginal beta
   hdis <- do.call(cbind, lapply(seq_along(alpha_prime), function(i) {
     hdi <- qbeta(c((1 - cred.int.level) / 2, (1 - ((1 - cred.int.level) / 2))),
       alpha_prime[i], sum(alpha_prime[-i]))
     return(hdi)
   }))
 
-  #' `calculate HDE`
-  #' note that this might change some since we could have HDE/HDI for every
-  #' marginal beta
+  #* `calculate HDE`
+  #* note that this might change some since we could have HDE/HDI for every
+  #* marginal beta
   hdes <- do.call(cbind, lapply(seq_along(alpha_prime), function(i) {
     hde <- .betaHDE(alpha_prime[i], sum(alpha_prime[-i]))
     return(hde)
   }))
   colnames(hdes) <- colnames(hdis) <- colnames(marginal_dens) <- colnames(marginal_post) <- nms
 
-  #' `Store summary`
-  #' note that this might change some since we could have HDE/HDI for every
-  #' marginal beta
+  #* `Store summary`
+  #* note that this might change some since we could have HDE/HDI for every
+  #* marginal beta
   out <- list()
   out$summary <- data.frame(
     HDE_1 = as.numeric(hdes), HDI_1_low = hdis[1, ], HDI_1_high = hdis[2, ]
@@ -125,10 +99,10 @@
   out$posteriorDraws <- marginal_post
   out$pdf <- marginal_dens
 
-  #' `Save s1 data for plotting`
-  #' note that this might change significantly from other distributions
-  #' because we want to test/plot the marginal betas to keep the dimensions
-  #' low enough to understand (and make the math tractable?)
+  #* `Save s1 data for plotting`
+  #* note that this might change significantly from other distributions
+  #* because we want to test/plot the marginal betas to keep the dimensions
+  #* low enough to understand (and make the math tractable?)
   out$plot_list <- list(
     "range" = range(support),
     "ddist_fun" = "stats::dbeta",
@@ -163,28 +137,9 @@
 }
 
 
-#' Might define the plot helper functions here as well, I don't think I have anything else
-#' that is quite in this situation for how it's using marginals, that's why I have the
-#' bivariate options for some distributions but a shared bivariate plotting function.
-#' Overall example would be drawn from `conjugate_plot::.conj_general_plot` and might just
-#' be pre-processing.
-#' Still need to decide how those things are done, thinking that it would be specified in
-#' the hypothesis, so instead of `(s1) ">" (s2)` I'd have the input be `(s1)K_1 > (s2)K_3`?
-#' I think that if there are 2 samples then it would assume your hypothesis has to be between
-#' the samples (previous example) but if there is only 1 sample that you're comparing categories
-#' in that sample.
-#' 
-#' This will take some real review of how hypotheses are handled to implement.
-#' Currently the hypothesis is passed to `.post.prob.from.pdfs` and is just used for logic in picking
-#' which comparison to make between PDFs, so I could have the dirichlet helper parse the hypothesis
-#' string and only return PDFs of the marginal betas that are requested? I think then I could return
-#' two for each sample based on the hypothesis, then just always take the first from s1 and the
-#' second from s2? That would mean I'd need a helper higher upstream to parse that hypothesis since
-#' by the time the `.conj_distHelper` function sees data it only has 1 sample.
-#' Not 100% sure that works right now but that feels reasonable.
-
 #' @keywords internal
 #' @noRd
+
 .multinomial.parse.hypothesis <- function(hypothesis) {
   g1 <- trimws(gsub("(.*?)([>|<|=|!]{1,2})(.*)", "\\1", hypothesis))
   g2 <- trimws(gsub("(.*?)([>|<|=|!]{1,2})(.*)", "\\3", hypothesis))
@@ -198,6 +153,7 @@
 
 #' @keywords internal
 #' @noRd
+
 .multinomial.pdf.handling <- function(sample_results, hypothesis) {
   parsed <- .multinomial.parse.hypothesis(hypothesis)
   g1 <- parsed[[1]]
@@ -228,11 +184,6 @@
 #' @keywords internal
 #' @noRd
 .multinomial.conj.plot.format <- function(res) {
-#' not immediately sure how this should work yet with the hypothesis being shuttled around
-  #' basically what I want to have happen is that this selects only the PDF/posteriors from the
-  #' groups that are in the hypothesis and plots those as betas. How to pass the hypothesis
-  #' information is tricky. Could label it earlier on? Could have an extra attribute in `res`
-  #' to store the formatted hypothesis?
   hypl <- .multinomial.parse.hypothesis(res$call$hypothesis)
   g1 <- hypl[[1]]
   g2 <- hypl[[2]]
@@ -285,7 +236,7 @@
   g1 <- hypl[[1]]
   g2 <- hypl[[2]]
   post1 <- sample_results[[1]]$posteriorDraws[, g1]
-  #' for ROPE comparison everything expects 2L results, so force that.
+  #* for ROPE comparison everything expects 2L results, so force that.
   if (length(sample_results) == 1) {
     post2 <- sample_results[[1]]$posteriorDraws[, g2]
     new_sample_results <- list(sample_results, sample_results)
@@ -297,15 +248,3 @@
   new_sample_results[[2]]$posteriorDraws <- post2
   return(new_sample_results)
 }
-
-#' @param s_res results from conjugate function thus far, currently the plot_list
-#' (for the distribution function name and values) element is all that is used.
-#' Internally this object is called `sample_results` in conjugate and only has
-#' one sample at a time passed to this function. Here it just needs to be
-#' formatted to only have the relevant pieces.
-#' @keywords internal
-#' @noRd
-
-.multinomial.bayes.factor.format <- function(s_res) {
-
-}
diff --git a/man/conjugate.Rd b/man/conjugate.Rd
index 152fc8ab..cd54c9e6 100644
--- a/man/conjugate.Rd
+++ b/man/conjugate.Rd
@@ -9,15 +9,14 @@ conjugate(
   s2 = NULL,
   method = c("t", "gaussian", "beta", "binomial", "lognormal", "lognormal2", "poisson",
     "negbin", "vonmises", "vonmises2", "uniform", "pareto", "gamma", "bernoulli",
-    "exponential", "bivariate_uniform", "bivariate_gaussian", "bivariate_lognormal"),
+    "exponential", "bivariate_uniform", "bivariate_gaussian", "bivariate_lognormal",
+    "multinomial"),
   priors = NULL,
-  plot = NULL,
   rope_range = NULL,
   rope_ci = 0.89,
   cred.int.level = 0.89,
   hypothesis = "equal",
-  bayes_factor = NULL,
-  support = NULL
+  bayes_factor = NULL
 )
 }
 \arguments{
@@ -34,9 +33,9 @@ This sample is shown in red if plotted.}
 This sample is shown in blue if plotted.}
 
 \item{method}{The distribution/method to use.
-Currently "t", "gaussian", "beta", "binomial", "lognormal", "lognormal2", "poisson",
-"negbin" (negative binomial), "uniform", "pareto", "gamma", "bernoulli", "exponential",
-"vonmises", and "vonmises2" are supported.
+Currently "bernoulli", "beta", "binomial", "exponential", "gamma", "gaussian",
+"lognormal", "lognormal2", "multinomial", "negbin" (negative binomial), "pareto",
+"poisson", "t", "uniform", "vonmises", and "vonmises2" are supported.
 The count (binomial, poisson and negative binomial), bernoulli, exponential,
 and pareto distributions are only implemented for single value traits due to their updating
 and/or the nature of the input data.
@@ -58,8 +57,6 @@ These names vary by method (see details).
  The \code{posterior} part of output can also be recycled as a new prior if Bayesian
  updating is appropriate for your use.}
 
-\item{plot}{deprecated, use \code{plot} method instead.}
-
 \item{rope_range}{Optional vector specifying a region of practical equivalence.
 This interval is considered practically equivalent to no effect.
 Kruschke (2018) suggests c(-0.1, 0.1) as a broadly reasonable ROPE for standardized parameters.
@@ -76,14 +73,22 @@ and \doi{10.1037/a0029146}.}
 in computing HDI for samples, defaults to 0.89.}
 
 \item{hypothesis}{Direction of a hypothesis if two samples are provided.
-Options are "unequal", "equal", "greater", and "lesser",
- read as "sample1 greater than sample2".}
+ Options are "unequal", "equal", "greater", and "lesser",
+  read as "sample1 greater than sample2". For the "multinomial" method the hypothesis
+should be specified as "Group1 >|<|==|!= Group2" and comparisons will be made using the marginal
+Beta distributions. If s2 is supplied then the hypothesis is read as
+"Group1 (from s1) >|<|==|!= Group2 (from s2)", if s2 is not supplied then both groups are taken
+from s1. For the multinomial method groups should be specified, so the hypothesis is written as
+"group1 > group2", where "group1" and "group2" would be informed by s1 unless s2 is provided in
+which case "group1" would be from s1 and "group2" would be from s2.}
 
 \item{bayes_factor}{Optional point or interval to evaluate bayes factors on. Note that this
 generally only makes sense to use if you have informative priors where the change in odds between
 prior and posterior is meaningful about the data. If this is non-NULL then columns of bayes factors
 are added to the summary output. Note these are only implemented for univariate distributions.}
 
+\item{plot}{deprecated, use \code{plot} method instead.}
+
 \item{support}{Deprecated}
 }
 \value{
@@ -173,6 +178,12 @@ Prior distributions default to be weakly informative and in some cases you may w
     kappa from the data and updates the kappa prior as a weighted average of the data and the prior.
     The mu parameter is then updated per Von-Mises conjugacy.
     }
+    \item{\strong{"multinomial": }
+    \code{list(alpha = list("alpha" = rep(1/N_groups, N_groups))}, where alpha is the concentration
+    vector of the conjugate dirichlet distribution. For the multinomial method hypotheses are
+    specified with group names, so instead of "equal" the hypothesis could be
+    "genotypeX == genotypeY".
+    }
     \item{\strong{"bivariate_uniform": }
     \code{list(location_l = 1, location_u = 2, scale = 1)}, where scale is the
     shared scale parameter of the pareto distributed upper and lower boundaries and location l and u
@@ -302,6 +313,26 @@ negbin_sv_ex <- conjugate(
   cred.int.level = 0.89, hypothesis = "equal"
 )
 
+# dirichlet-multinomial sv example
+
+dm_sv_ex <- conjugate(
+  s1 = list("A" = 10, "B" = 10, "C" = 5),
+  s2 = list("A" = 4, "B" = 12, "C" = 9),
+  method = "multinomial",
+  hypothesis = "A > A",
+  rope_range = c(-0.1, 0.1)
+)
+
+# dirichlet-multinomial mv example
+
+dm_mv_ex <- conjugate(
+  s1 = data.frame("A" = c(5,5), "B" = c(5, 5), "C" = c(2,3)),
+  s2 = data.frame("A" = c(2,2), "B" = c(7, 5), "C" = c(8,1)),
+  method = "multinomial",
+  hypothesis = "A > A",
+  rope_range = c(-0.1, 0.1)
+)
+
 # von mises mv example
 
 mv_gauss <- mvSim(

From 530cdc25357d11e418020db4a7406b97554a201e Mon Sep 17 00:00:00 2001
From: Josh Sumner <51797700+joshqsumner@users.noreply.github.com>
Date: Thu, 23 Apr 2026 13:50:32 -0500
Subject: [PATCH 19/25] linting

---
 R/conjugate.R                    |  6 ++++--
 R/conjugate_class.R              | 16 ++++++++++------
 R/conjugate_multinomialHelpers.R | 23 ++++++++++++-----------
 3 files changed, 26 insertions(+), 19 deletions(-)

diff --git a/R/conjugate.R b/R/conjugate.R
index 7dfc87ba..b6524373 100644
--- a/R/conjugate.R
+++ b/R/conjugate.R
@@ -408,8 +408,10 @@ conjugate <- function(s1 = NULL, s2 = NULL,
     mult_prob <- .multinomial.pdf.handling(sample_results, hypothesis)
     out$summary <- cbind(
       out$summary,
-      data.frame("hyp" = mult_prob$hyp,
-        "post.prob" = as.numeric(mult_prob$pdf.handling.output$post.prob))
+      data.frame(
+        "hyp" = mult_prob$hyp,
+        "post.prob" = as.numeric(mult_prob$pdf.handling.output$post.prob)
+      )
     )
   } else if (!is.null(s2) && method[1] != "multinomial") {
     postProbRes <- .pdf.handling(sample_results[[1]]$pdf, sample_results[[2]]$pdf, hypothesis)
diff --git a/R/conjugate_class.R b/R/conjugate_class.R
index 53e6d7a9..4f831984 100644
--- a/R/conjugate_class.R
+++ b/R/conjugate_class.R
@@ -112,10 +112,12 @@ print.conjugatesummary <- function(x, ...) {
       prior_statement <- paste0(
         "Sample ", i, " Prior ", method_list[[1]], "(",
         paste(
-          paste(names(prior),    
+          paste(
+            names(prior),
             unlist(lapply(names(prior), function(nm) {
-              paste(round(prior[[nm]], 3), collapse = ", ")
-            })), sep = " = "),
+              return(paste(round(prior[[nm]], 3), collapse = ", "))
+            })), sep = " = "
+          ),
           collapse = ", "
         ),
         ")\n"
@@ -123,10 +125,12 @@ print.conjugatesummary <- function(x, ...) {
       posterior_statement <- paste0(
         "\tPosterior ", method_list[[1]], "(",
         paste(
-          paste(names(post),    
+          paste(
+            names(post),
             unlist(lapply(names(post), function(nm) {
-              paste(round(post[[nm]], 3), collapse = ", ")
-            })), sep = " = "),
+              return(paste(round(post[[nm]], 3), collapse = ", "))
+            })), sep = " = "
+          ),
           collapse = ", "
         ),
         ")\n"
diff --git a/R/conjugate_multinomialHelpers.R b/R/conjugate_multinomialHelpers.R
index dc6a5c09..5811ba4e 100644
--- a/R/conjugate_multinomialHelpers.R
+++ b/R/conjugate_multinomialHelpers.R
@@ -5,7 +5,9 @@
 #' Note that this function returns posterior draws and densities from the marginal beta
 #' distributions, NOT from the dirichlet distribution.
 #'
-#' @param s1 A named vector/list of numerics drawn from a multinomial distribution. A data.frame will trigger the mv option, which will take column sums of the data then pass to the sv method.
+#' @param s1 A named vector/list of numerics drawn from a multinomial distribution.
+#' A data.frame will trigger the mv option,
+#' which will take column sums of the data then pass to the sv method.
 #' @examples
 #' out <- .conj_dirichlet_sv(
 #'   s1 = list("A" = 10, "B" = 10, "C" = 5),
@@ -22,7 +24,7 @@
 #' hypothesis = "A > B"
 #' sample_number = 1
 #'
-#' 
+#'
 #' @details
 #' See Examples 1.4, 1.6, and 1.7 for thoughts on default dirichlet prior here
 #' https://arxiv.org/pdf/1504.02689 , updating rule defined in
@@ -32,8 +34,8 @@
 #' @noRd
 
 .conj_multinomial_sv <- function(s1 = NULL, priors = NULL,
-                               support = NULL, cred.int.level = NULL,
-                               calculatingSupport = FALSE, ...) {
+                                 support = NULL, cred.int.level = NULL,
+                                 calculatingSupport = FALSE, ...) {
   nms <- names(s1)
   s1 <- unlist(s1)
 
@@ -42,8 +44,6 @@
   if (is.null(support) && calculatingSupport) {
     return(c(0.0001, 0.9999))
   }
-  #* `Check stopping conditions (similar to binomial)`
-  
   #* `make default prior if none provided`
   if (is.null(priors)) {
     priors <- list("alpha" = rep(1 / length(s1), length(s1)))
@@ -72,8 +72,10 @@
   #* note that this might change some since we could have HDE/HDI for every
   #* marginal beta
   hdis <- do.call(cbind, lapply(seq_along(alpha_prime), function(i) {
-    hdi <- qbeta(c((1 - cred.int.level) / 2, (1 - ((1 - cred.int.level) / 2))),
-      alpha_prime[i], sum(alpha_prime[-i]))
+    hdi <- qbeta(
+      c((1 - cred.int.level) / 2, (1 - ((1 - cred.int.level) / 2))),
+      alpha_prime[i], sum(alpha_prime[-i])
+    )
     return(hdi)
   }))
 
@@ -127,8 +129,8 @@
 #' @noRd
 
 .conj_multinomial_mv <- function(s1 = NULL, priors = NULL,
-                               support = NULL, cred.int.level = NULL,
-                               calculatingSupport = FALSE, ...) {
+                                 support = NULL, cred.int.level = NULL,
+                                 calculatingSupport = FALSE, ...) {
   nms <- colnames(s1)
   s1 <- colSums(s1)
   names(s1) <- nms
@@ -187,7 +189,6 @@
   hypl <- .multinomial.parse.hypothesis(res$call$hypothesis)
   g1 <- hypl[[1]]
   g2 <- hypl[[2]]
-  hyp <- hypl[[3]]
   new_res <- res
   if (length(res$plot_parameters) == 2) {
     # keep only relevant marginal beta parameters

From 708cf27b2a8fea3562e2ff8991d8e88b6cba2b76 Mon Sep 17 00:00:00 2001
From: Josh Sumner <51797700+joshqsumner@users.noreply.github.com>
Date: Thu, 23 Apr 2026 14:17:15 -0500
Subject: [PATCH 20/25] removing support and plot from examples

---
 R/conjugate.R                      | 12 +++++-------
 man/conjugate.Rd                   | 12 ++++--------
 tests/testthat/test-sv-conjugate.R |  2 +-
 3 files changed, 10 insertions(+), 16 deletions(-)

diff --git a/R/conjugate.R b/R/conjugate.R
index b6524373..c263d312 100644
--- a/R/conjugate.R
+++ b/R/conjugate.R
@@ -37,7 +37,6 @@
 #'  By default this is NULL and weak priors (generally jeffrey's priors) are used.
 #'  The \code{posterior} part of output can also be recycled as a new prior if Bayesian
 #'  updating is appropriate for your use.
-#' @param plot deprecated, use \code{plot} method instead.
 #' @param rope_range Optional vector specifying a region of practical equivalence.
 #' This interval is considered practically equivalent to no effect.
 #' Kruschke (2018) suggests c(-0.1, 0.1) as a broadly reasonable ROPE for standardized parameters.
@@ -63,7 +62,6 @@
 #' generally only makes sense to use if you have informative priors where the change in odds between
 #' prior and posterior is meaningful about the data. If this is non-NULL then columns of bayes factors
 #' are added to the summary output. Note these are only implemented for univariate distributions.
-#' @param support Deprecated
 #'
 #' @import bayestestR
 #' @import extraDistr
@@ -169,7 +167,7 @@
 #'   s1 = mv_ln[1:30, -1], s2 = mv_ln[31:60, -1], method = "lognormal",
 #'   priors = list(mu = 5, sd = 2),
 #'   rope_range = c(-40, 40), rope_ci = 0.89,
-#'   cred.int.level = 0.89, hypothesis = "equal", support = NULL
+#'   cred.int.level = 0.89, hypothesis = "equal"
 #' )
 #'
 #' # lognormal sv
@@ -178,7 +176,7 @@
 #'   method = "lognormal",
 #'   priors = list(mu = 5, sd = 2),
 #'   rope_range = NULL, rope_ci = 0.89,
-#'   cred.int.level = 0.89, hypothesis = "equal", support = NULL
+#'   cred.int.level = 0.89, hypothesis = "equal"
 #' )
 #'
 #' # Z test mv example
@@ -195,7 +193,7 @@
 #'   s1 = mv_gauss[1:30, -1], s2 = mv_gauss[31:60, -1], method = "gaussian",
 #'   priors = list(mu = 30, sd = 10),
 #'   rope_range = c(-25, 25), rope_ci = 0.89,
-#'   cred.int.level = 0.89, hypothesis = "equal", support = NULL
+#'   cred.int.level = 0.89, hypothesis = "equal"
 #' )
 #'
 #' # T test sv example with two different priors
@@ -204,7 +202,7 @@
 #'   s1 = rnorm(10, 50, 10), s2 = rnorm(10, 60, 12), method = "t",
 #'   priors = list(list(mu = 40, sd = 10), list(mu = 45, sd = 8)),
 #'   rope_range = c(-5, 8), rope_ci = 0.89,
-#'   cred.int.level = 0.89, hypothesis = "equal", support = NULL
+#'   cred.int.level = 0.89, hypothesis = "equal"
 #' )
 #'
 #' # beta mv example
@@ -668,7 +666,7 @@ conjugate <- function(s1 = NULL, s2 = NULL,
 #' @keywords internal
 #' @noRd
 
-.pdf.handling <- function(method, pdf1, pdf2, hypothesis) {
+.pdf.handling <- function(pdf1, pdf2, hypothesis) {
   if (is.list(pdf1) && is.list(pdf2)) {
     pdf.handling.output <- as.data.frame(do.call(rbind, lapply(seq_along(pdf1), function(i) {
       pdf <- .post.prob.from.pdfs(pdf1[[i]], pdf2[[i]], hypothesis)
diff --git a/man/conjugate.Rd b/man/conjugate.Rd
index cd54c9e6..897c7d96 100644
--- a/man/conjugate.Rd
+++ b/man/conjugate.Rd
@@ -86,10 +86,6 @@ which case "group1" would be from s1 and "group2" would be from s2.}
 generally only makes sense to use if you have informative priors where the change in odds between
 prior and posterior is meaningful about the data. If this is non-NULL then columns of bayes factors
 are added to the summary output. Note these are only implemented for univariate distributions.}
-
-\item{plot}{deprecated, use \code{plot} method instead.}
-
-\item{support}{Deprecated}
 }
 \value{
 A \link{conjugate-class} object with slots including:
@@ -214,7 +210,7 @@ ln_mv_ex <- conjugate(
   s1 = mv_ln[1:30, -1], s2 = mv_ln[31:60, -1], method = "lognormal",
   priors = list(mu = 5, sd = 2),
   rope_range = c(-40, 40), rope_ci = 0.89,
-  cred.int.level = 0.89, hypothesis = "equal", support = NULL
+  cred.int.level = 0.89, hypothesis = "equal"
 )
 
 # lognormal sv
@@ -223,7 +219,7 @@ ln_sv_ex <- conjugate(
   method = "lognormal",
   priors = list(mu = 5, sd = 2),
   rope_range = NULL, rope_ci = 0.89,
-  cred.int.level = 0.89, hypothesis = "equal", support = NULL
+  cred.int.level = 0.89, hypothesis = "equal"
 )
 
 # Z test mv example
@@ -240,7 +236,7 @@ gauss_mv_ex <- conjugate(
   s1 = mv_gauss[1:30, -1], s2 = mv_gauss[31:60, -1], method = "gaussian",
   priors = list(mu = 30, sd = 10),
   rope_range = c(-25, 25), rope_ci = 0.89,
-  cred.int.level = 0.89, hypothesis = "equal", support = NULL
+  cred.int.level = 0.89, hypothesis = "equal"
 )
 
 # T test sv example with two different priors
@@ -249,7 +245,7 @@ gaussianMeans_sv_ex <- conjugate(
   s1 = rnorm(10, 50, 10), s2 = rnorm(10, 60, 12), method = "t",
   priors = list(list(mu = 40, sd = 10), list(mu = 45, sd = 8)),
   rope_range = c(-5, 8), rope_ci = 0.89,
-  cred.int.level = 0.89, hypothesis = "equal", support = NULL
+  cred.int.level = 0.89, hypothesis = "equal"
 )
 
 # beta mv example
diff --git a/tests/testthat/test-sv-conjugate.R b/tests/testthat/test-sv-conjugate.R
index be127b45..1c837059 100644
--- a/tests/testthat/test-sv-conjugate.R
+++ b/tests/testthat/test-sv-conjugate.R
@@ -98,7 +98,7 @@ test_that("conjugate single value lognormal works", {
   expect_warning(
     out <- conjugate(
       s1 = s1, s2 = s2,
-      method = "lognormal", priors = NULL, plot = TRUE,
+      method = "lognormal", priors = NULL,
       rope_range = c(-1, 1), rope_ci = 0.89, cred.int.level = 0.89,
       hypothesis = "equal", bayes_factor = 5
     )

From 3081fd1779d2a4fd8a9fbbdcdec647015078d0d2 Mon Sep 17 00:00:00 2001
From: Josh Sumner <51797700+joshqsumner@users.noreply.github.com>
Date: Thu, 23 Apr 2026 14:20:12 -0500
Subject: [PATCH 21/25] remove plot arg from barg example

---
 R/barg.R    | 2 +-
 man/barg.Rd | 2 +-
 2 files changed, 2 insertions(+), 2 deletions(-)

diff --git a/R/barg.R b/R/barg.R
index 174d7e9a..d27dba2e 100644
--- a/R/barg.R
+++ b/R/barg.R
@@ -129,7 +129,7 @@
 #'     list(mu = 10, sd = 2),
 #'     list(mu = 10, sd = 2)
 #'   ),
-#'   plot = FALSE, rope_range = c(-8, 8), rope_ci = 0.89,
+#'   rope_range = c(-8, 8), rope_ci = 0.89,
 #'   cred.int.level = 0.89, hypothesis = "unequal",
 #'   bayes_factor = c(50, 55)
 #' )
diff --git a/man/barg.Rd b/man/barg.Rd
index 941be0cf..ccacd7d4 100644
--- a/man/barg.Rd
+++ b/man/barg.Rd
@@ -136,7 +136,7 @@ x <- conjugate(
     list(mu = 10, sd = 2),
     list(mu = 10, sd = 2)
   ),
-  plot = FALSE, rope_range = c(-8, 8), rope_ci = 0.89,
+  rope_range = c(-8, 8), rope_ci = 0.89,
   cred.int.level = 0.89, hypothesis = "unequal",
   bayes_factor = c(50, 55)
 )

From 3e60a8b8a7828d31e95a85739332712d1e3fd4f0 Mon Sep 17 00:00:00 2001
From: Josh Sumner <51797700+joshqsumner@users.noreply.github.com>
Date: Fri, 24 Apr 2026 11:56:44 -0500
Subject: [PATCH 22/25] deepsource

---
 R/barg.R                           |  2 +-
 R/stat_brms_model.R                | 12 ++++++++----
 R/stat_nlme_model.R                |  2 +-
 man/barg.Rd                        |  2 +-
 tests/testthat/test-brmsModels.R   |  2 --
 tests/testthat/test-growthModels.R |  4 ++--
 tests/testthat/test-sv-conjugate.R | 12 +++++-------
 7 files changed, 18 insertions(+), 18 deletions(-)

diff --git a/R/barg.R b/R/barg.R
index d27dba2e..a40d766a 100644
--- a/R/barg.R
+++ b/R/barg.R
@@ -118,7 +118,7 @@
 #'   iter = 600, cores = 1, chains = 1, backend = "cmdstanr",
 #'   sample_prior = "only" # only sampling from prior for speed
 #' )
-#' barg(fit_test, ss)
+#' b <- barg(fit_test, ss)
 #' fit_2 <- fit_test
 #' fit_list <- list(fit_test, fit_2)
 #' x <- barg(fit_list, list(ss, ss))
diff --git a/R/stat_brms_model.R b/R/stat_brms_model.R
index 7a9b6122..ce8a8d5a 100644
--- a/R/stat_brms_model.R
+++ b/R/stat_brms_model.R
@@ -178,7 +178,7 @@ statBrmsMod <- ggplot2::ggproto("StatBrm", Stat,
 
 statBrmsStaticMod <- ggplot2::ggproto("StatStaticBrm", Stat,
   # `specify that there will be extra params`
-  extra_params = c("na.rm", "fit", "parsed_form", "probs"),
+  extra_params = c("na.rm", "fit", "parsed_form", "probs", "hierarchy_value"),
   # `data setup function`
   setup_data = function(data, params) {
     if ("data" %in% names(params$parsed_form)) {
@@ -190,6 +190,8 @@ statBrmsStaticMod <- ggplot2::ggproto("StatStaticBrm", Stat,
       data <- plyr::join(mod_data, data, type = "left", match = "all", by = "x")
       if (length(unique(data$PANEL)) > 1) {
         data <- data[data$PANEL == as.numeric(as.factor(data$MOD_GROUP)), ]
+      } else {
+        data$PANEL <- 1
       }
     }
     return(data)
@@ -268,8 +270,9 @@ statBrmsStaticMod <- ggplot2::ggproto("StatStaticBrm", Stat,
     longPreds$xmax <- longPreds$numericGroup + c(0.45 * (1 - longPreds$q))
 
     # select columns and rename
-    grpdf <- longPreds[, c(group, "Estimate", "PANEL", "q", "ymin", "ymax", "xmin", "xmax")]
-    colnames(grpdf) <- c("MOD_GROUP", "y", "PANEL", "Cred.Int", "ymin", "ymax", "xmin", "xmax")
+    grpdf <- longPreds[, c(group,
+      "Estimate", "PANEL", "q", "ymin", "ymax", "xmin", "xmax", "numericGroup")]
+    colnames(grpdf) <- c("MOD_GROUP", "y", "PANEL", "Cred.Int", "ymin", "ymax", "xmin", "xmax", "grp")
     return(grpdf)
   },
   # set defaults for several aesthetics, all have to be after stat is calculated
@@ -280,6 +283,7 @@ statBrmsStaticMod <- ggplot2::ggproto("StatStaticBrm", Stat,
     xmax = after_stat(xmax),
     fill = after_stat(Cred.Int),
     alpha = after_stat(0.5),
-    group = after_stat(x)
+    group = after_stat(grp),
+    x = after_stat(MOD_GROUP)
   )
 )
diff --git a/R/stat_nlme_model.R b/R/stat_nlme_model.R
index a84f941a..f1509715 100644
--- a/R/stat_nlme_model.R
+++ b/R/stat_nlme_model.R
@@ -51,7 +51,7 @@ statNlmeMod <- ggplot2::ggproto("StatNlme", Stat,
         nrow(mod_data), "\n\n", paste(
           do.call(
             cbind,
-            lapply(mod_data, \(x) paste(unique(x), collapse = ", "))
+            lapply(mod_data, function(x) paste(unique(x), collapse = ", "))
           ),
           collapse = " -- "
         )
diff --git a/man/barg.Rd b/man/barg.Rd
index ccacd7d4..ad453617 100644
--- a/man/barg.Rd
+++ b/man/barg.Rd
@@ -125,7 +125,7 @@ fit_test <- fitGrowth(ss,
   iter = 600, cores = 1, chains = 1, backend = "cmdstanr",
   sample_prior = "only" # only sampling from prior for speed
 )
-barg(fit_test, ss)
+b <- barg(fit_test, ss)
 fit_2 <- fit_test
 fit_list <- list(fit_test, fit_2)
 x <- barg(fit_list, list(ss, ss))
diff --git a/tests/testthat/test-brmsModels.R b/tests/testthat/test-brmsModels.R
index 62632e94..64e71227 100644
--- a/tests/testthat/test-brmsModels.R
+++ b/tests/testthat/test-brmsModels.R
@@ -21,13 +21,11 @@ test_that("Logistic brms model pipeline", {
     df = simdf, type = "brms"
   )
   expect_equal(ss$prior$nlpar, c("", "", "A", "B", "C"))
-
   fit <- fitGrowth(ss,
     backend = "cmdstanr", iter = 500, chains = 1, cores = 1,
     refresh = 0, silent = 2
   )
   expect_s3_class(fit, "brmsfit")
-
   plot <- growthPlot(fit = fit, form = ss$pcvrForm, df = ss$df)
   expect_s3_class(plot, "ggplot")
   p <- ggplot() + stat_growthss(fit  = fit, ss = ss)
diff --git a/tests/testthat/test-growthModels.R b/tests/testthat/test-growthModels.R
index 0c15adf7..2b8e4664 100644
--- a/tests/testthat/test-growthModels.R
+++ b/tests/testthat/test-growthModels.R
@@ -71,7 +71,7 @@ test_that("Test Logistic nls modeling", {
       3.34892124655795, 3.34892124655795
     )
   )
-  invisible(print(ss))
+  x <- invisible(print(ss))
   fit <- fitGrowth(ss)
   expect_s3_class(fit, "nls")
 
@@ -131,7 +131,7 @@ test_that("Test Logistic nlme modeling", {
       3.34892124655795, 3.34892124655795
     )
   )
-  invisible(print(ss))
+  x <- invisible(print(ss))
   fit <- suppressWarnings(fitGrowth(ss))
   expect_s3_class(fit, "nlme")
 
diff --git a/tests/testthat/test-sv-conjugate.R b/tests/testthat/test-sv-conjugate.R
index 1c837059..7ff8422e 100644
--- a/tests/testthat/test-sv-conjugate.R
+++ b/tests/testthat/test-sv-conjugate.R
@@ -95,13 +95,11 @@ test_that("conjugate single value lognormal works", {
   set.seed(123)
   s1 <- rlnorm(100, log(130), log(1.3))
   s2 <- rlnorm(100, log(100), log(2))
-  expect_warning(
-    out <- conjugate(
-      s1 = s1, s2 = s2,
-      method = "lognormal", priors = NULL,
-      rope_range = c(-1, 1), rope_ci = 0.89, cred.int.level = 0.89,
-      hypothesis = "equal", bayes_factor = 5
-    )
+  out <- conjugate(
+    s1 = s1, s2 = s2,
+    method = "lognormal", priors = NULL,
+    rope_range = c(-1, 1), rope_ci = 0.89, cred.int.level = 0.89,
+    hypothesis = "equal", bayes_factor = 5
   )
   expect_equal(out$summary$post.prob, 0.5980101, tolerance = 1e-6)
   expect_equal(out$summary$rope_prob, 0.1113358, tolerance = 1e-6)

From 6309fca5eba84f42e0ccb4b7a4d80b6ea51d5d1c Mon Sep 17 00:00:00 2001
From: Josh Sumner <51797700+joshqsumner@users.noreply.github.com>
Date: Fri, 24 Apr 2026 14:16:33 -0500
Subject: [PATCH 23/25] lint

---
 R/stat_brms_model.R | 6 ++----
 1 file changed, 2 insertions(+), 4 deletions(-)

diff --git a/R/stat_brms_model.R b/R/stat_brms_model.R
index ce8a8d5a..7be6819b 100644
--- a/R/stat_brms_model.R
+++ b/R/stat_brms_model.R
@@ -270,9 +270,8 @@ statBrmsStaticMod <- ggplot2::ggproto("StatStaticBrm", Stat,
     longPreds$xmax <- longPreds$numericGroup + c(0.45 * (1 - longPreds$q))
 
     # select columns and rename
-    grpdf <- longPreds[, c(group,
-      "Estimate", "PANEL", "q", "ymin", "ymax", "xmin", "xmax", "numericGroup")]
-    colnames(grpdf) <- c("MOD_GROUP", "y", "PANEL", "Cred.Int", "ymin", "ymax", "xmin", "xmax", "grp")
+    grpdf <- longPreds[, c(group, "Estimate", "PANEL", "q", "ymin", "ymax", "xmin", "xmax")]
+    colnames(grpdf) <- c("MOD_GROUP", "y", "PANEL", "Cred.Int", "ymin", "ymax", "xmin", "xmax")
     return(grpdf)
   },
   # set defaults for several aesthetics, all have to be after stat is calculated
@@ -283,7 +282,6 @@ statBrmsStaticMod <- ggplot2::ggproto("StatStaticBrm", Stat,
     xmax = after_stat(xmax),
     fill = after_stat(Cred.Int),
     alpha = after_stat(0.5),
-    group = after_stat(grp),
     x = after_stat(MOD_GROUP)
   )
 )

From 29aad02c110c4255b66ef9782f57beffd57ed526 Mon Sep 17 00:00:00 2001
From: Josh Sumner <51797700+joshqsumner@users.noreply.github.com>
Date: Fri, 24 Apr 2026 15:17:00 -0500
Subject: [PATCH 24/25] fixing 1 sample plot

---
 R/conjugate_class.R                |  2 +-
 R/conjugate_multinomialHelpers.R   | 18 ++++++++++-----
 tests/testthat/test-sv-conjugate.R | 36 ++++++++++++++++++++++++++++++
 3 files changed, 49 insertions(+), 7 deletions(-)

diff --git a/R/conjugate_class.R b/R/conjugate_class.R
index 4f831984..ecda19ac 100644
--- a/R/conjugate_class.R
+++ b/R/conjugate_class.R
@@ -171,7 +171,7 @@ print.conjugatesummary <- function(x, ...) {
       100 * rope_ci, "% Credible Interval is ",
       100 * round(x$summary$rope_prob[1], 5),
       "% with an average difference of ",
-      round(x$summary$HDE_rope, 3)
+      round(x$summary$HDE_rope[1], 3)
     )
     cat(rope_message)
     cat("\n\n")
diff --git a/R/conjugate_multinomialHelpers.R b/R/conjugate_multinomialHelpers.R
index 5811ba4e..fba2df1f 100644
--- a/R/conjugate_multinomialHelpers.R
+++ b/R/conjugate_multinomialHelpers.R
@@ -200,7 +200,7 @@
     new_res$plot_parameters[[1]]$priors$shape1 <- res$plot_parameters[[1]]$priors$shape1[g1]
     new_res$plot_parameters[[2]]$priors$shape1 <- res$plot_parameters[[2]]$priors$shape1[g2]
     new_res$plot_parameters[[1]]$priors$shape2 <- res$plot_parameters[[1]]$priors$shape2[g1]
-    new_res$plot_parameters[[2]]$priors$shape2 <- res$plot_parameters[[2]]$priors$shape2[g1]
+    new_res$plot_parameters[[2]]$priors$shape2 <- res$plot_parameters[[2]]$priors$shape2[g2]
     # subset summary to draw HDE/HDI correctly
     new_res$summary <- cbind(
       res$summary[g1, grepl("_1", colnames(res$summary))],
@@ -208,14 +208,20 @@
       res$summary[1, !grepl("[1|2]", colnames(res$summary))]
     )
   } else {
+    # duplicate plot parameters to 2L
+    new_res$plot_parameters <- list(new_res$plot_parameters[[1]], new_res$plot_parameters[[1]])
+    new_res$data <- list(new_res$data[[1]], new_res$data[[1]])
     # keep only relevant marginal beta parameters
-    new_res$plot_parameters[[1]]$parameters$shape1 <- res$plot_parameters[[1]]$parameters$shape1[g1]
-    new_res$plot_parameters[[1]]$parameters$shape1 <- res$plot_parameters[[1]]$parameters$shape2[g2]
+    new_res$plot_parameters[[1]]$parameters$shape1 <- new_res$plot_parameters[[1]]$parameters$shape1[g1]
+    new_res$plot_parameters[[2]]$parameters$shape1 <- new_res$plot_parameters[[2]]$parameters$shape1[g2]
+    new_res$plot_parameters[[1]]$parameters$shape2 <- new_res$plot_parameters[[1]]$parameters$shape2[g1]
+    new_res$plot_parameters[[2]]$parameters$shape2 <- new_res$plot_parameters[[2]]$parameters$shape2[g2]
     # keep only relevant marginal priors
-    new_res$plot_parameters[[1]]$priors$shape1 <- res$plot_parameters[[1]]$priors$shape1[g1]
-    new_res$plot_parameters[[1]]$priors$shape1 <- res$plot_parameters[[1]]$priors$shape2[g2]
+    new_res$plot_parameters[[1]]$priors$shape1 <- new_res$plot_parameters[[1]]$priors$shape1[g1]
+    new_res$plot_parameters[[2]]$priors$shape1 <- new_res$plot_parameters[[2]]$priors$shape1[g2]
+    new_res$plot_parameters[[1]]$priors$shape2 <- new_res$plot_parameters[[1]]$priors$shape2[g1]
+    new_res$plot_parameters[[2]]$priors$shape2 <- new_res$plot_parameters[[2]]$priors$shape2[g2]
     # subset summary to draw HDE/HDI correctly
-    ls()
     g1_cols <- res$summary[g1, grepl("_1", colnames(res$summary))]
     g2_cols <- res$summary[g2, grepl("_1", colnames(res$summary))]
     colnames(g2_cols) <- gsub("_1", "_2", colnames(g2_cols))
diff --git a/tests/testthat/test-sv-conjugate.R b/tests/testthat/test-sv-conjugate.R
index 7ff8422e..b4d91998 100644
--- a/tests/testthat/test-sv-conjugate.R
+++ b/tests/testthat/test-sv-conjugate.R
@@ -336,6 +336,42 @@ test_that("conjugate single value exponential works", {
   expect_s3_class(out, "conjugate")
 })
 
+test_that("conjugate single value Dirichlet-Multinomial works", {
+  set.seed(123)
+  s1 <- list(A = 10, B = 20, C = 30)
+  s2 <- list(A = 10, B = 20, C = 30)
+  out <- conjugate(
+    s1 = s1, s2 = s2, method = "multinomial",
+    priors = NULL,
+    rope_range = c(-0.15, 0.15), rope_ci = 0.89,
+    cred.int.level = 0.89, hypothesis = "A > B",
+    bayes_factor = 0.5
+  )
+  expect_equal(out$summary$post.prob[1], 0.9306264, tolerance = 1e-6)
+  expect_equal(out$summary$rope_prob[1], 0.411077407032918, tolerance = 1e-6)
+  expect_s3_class(out, "conjugate")
+  p <- plot(out)
+  expect_s3_class(p, "ggplot")
+})
+
+
+test_that("conjugate one sample single value Dirichlet-Multinomial works", {
+  set.seed(123)
+  s1 <- list(A = 10, B = 20, C = 30)
+  out <- conjugate(
+    s1 = s1, method = "multinomial",
+    priors = NULL,
+    rope_range = c(-0.15, 0.15), rope_ci = 0.89,
+    cred.int.level = 0.89, hypothesis = "A > B",
+    bayes_factor = 0.5
+  )
+  expect_equal(out$summary$post.prob[1], 0.9306264, tolerance = 1e-6)
+  expect_equal(out$summary$rope_prob[1], 0.4209639, tolerance = 1e-6)
+  expect_s3_class(out, "conjugate")
+  p <- plot(out)
+  expect_s3_class(p, "ggplot")
+})
+
 test_that("conjugate single value lognormal vs gaussian", {
   set.seed(123)
   s1 <- rlnorm(100, log(70), log(2))

From 9c302943a2e0ed01ea112e24c5aa213ca91706d7 Mon Sep 17 00:00:00 2001
From: Josh Sumner <51797700+joshqsumner@users.noreply.github.com>
Date: Mon, 27 Apr 2026 09:38:20 -0500
Subject: [PATCH 25/25] fix test

---
 tests/testthat/test-sv-conjugate.R | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/tests/testthat/test-sv-conjugate.R b/tests/testthat/test-sv-conjugate.R
index b4d91998..85a13fa6 100644
--- a/tests/testthat/test-sv-conjugate.R
+++ b/tests/testthat/test-sv-conjugate.R
@@ -348,7 +348,7 @@ test_that("conjugate single value Dirichlet-Multinomial works", {
     bayes_factor = 0.5
   )
   expect_equal(out$summary$post.prob[1], 0.9306264, tolerance = 1e-6)
-  expect_equal(out$summary$rope_prob[1], 0.411077407032918, tolerance = 1e-6)
+  expect_equal(out$summary$rope_prob[1], 0.4252331, tolerance = 1e-6)
   expect_s3_class(out, "conjugate")
   p <- plot(out)
   expect_s3_class(p, "ggplot")