diff --git a/R/barg.R b/R/barg.R
index 174d7e9a..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))
@@ -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/R/conjugate.R b/R/conjugate.R
index ac45eff4..c263d312 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.
@@ -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.
@@ -52,12 +51,17 @@
 #' 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. 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
 #' are added to the summary output. Note these are only implemented for univariate distributions.
-#' @param support Deprecated
 #'
 #' @import bayestestR
 #' @import extraDistr
@@ -127,6 +131,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
@@ -157,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
@@ -166,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
@@ -183,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
@@ -192,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
@@ -256,6 +266,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(
@@ -313,11 +343,11 @@ 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, 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 +363,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) {
@@ -357,7 +380,8 @@ 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)
@@ -368,7 +392,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) {
@@ -377,7 +401,17 @@ conjugate <- function(s1 = NULL, s2 = NULL,
     }
     return(s)
   }))
-  if (!is.null(s2)) {
+  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(
+      out$summary,
+      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)
     out$summary <- cbind(
       out$summary,
@@ -386,7 +420,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
   }
@@ -536,8 +570,9 @@ 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) {
+                       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))) {
@@ -547,6 +582,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_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_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)
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_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)
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_class.R b/R/conjugate_class.R
index 13faa9a8..ecda19ac 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
@@ -111,7 +112,12 @@ 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) {
+              return(paste(round(prior[[nm]], 3), collapse = ", "))
+            })), sep = " = "
+          ),
           collapse = ", "
         ),
         ")\n"
@@ -119,7 +125,12 @@ 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) {
+              return(paste(round(post[[nm]], 3), collapse = ", "))
+            })), sep = " = "
+          ),
           collapse = ", "
         ),
         ")\n"
@@ -160,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_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_multinomialHelpers.R b/R/conjugate_multinomialHelpers.R
new file mode 100644
index 00000000..fba2df1f
--- /dev/null
+++ b/R/conjugate_multinomialHelpers.R
@@ -0,0 +1,257 @@
+#' @description
+#' Internal function for calculating the dirichlet distribution underlying
+#' multinomially distributed single value traits (counts).
+#'
+#' 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),
+#'   cred.int.level = 0.95,
+#'   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
+#' 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,
+                                 support = NULL, cred.int.level = NULL,
+                                 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))
+  }
+  #* `make default prior if none provided`
+  if (is.null(priors)) {
+    priors <- list("alpha" = rep(1 / length(s1), length(s1)))
+  }
+  names(priors$alpha) <- nms
+
+  #* `update dirichlet prior with sufficient statistics (counts)`
+  alpha_prime <- priors$alpha + s1
+  #* `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 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
+  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
+  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" = setNames(unlist(lapply(seq_along(priors$alpha), function(i) {
+        return(sum(priors$alpha[-i]))
+      })), nms)
+    ),
+    "parameters" = list("shape1" = alpha_prime,
+      "shape2" = setNames(unlist(lapply(seq_along(alpha_prime), function(i) {
+        return(sum(alpha_prime[-i]))
+      })), nms)
+    ),
+    "given" = NULL
+  )
+  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)
+}
+
+
+#' @keywords internal
+#' @noRd
+
+.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],
+      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
+    )
+  )
+}
+
+
+#' @keywords internal
+#' @noRd
+.multinomial.conj.plot.format <- function(res) {
+  hypl <- .multinomial.parse.hypothesis(res$call$hypothesis)
+  g1 <- hypl[[1]]
+  g2 <- hypl[[2]]
+  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[g2]
+    # 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 {
+    # 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 <- 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 <- 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
+    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)
+}
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_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)
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`
diff --git a/R/stat_brms_model.R b/R/stat_brms_model.R
index 7a9b6122..7be6819b 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)
@@ -280,6 +282,6 @@ statBrmsStaticMod <- ggplot2::ggproto("StatStaticBrm", Stat,
     xmax = after_stat(xmax),
     fill = after_stat(Cred.Int),
     alpha = after_stat(0.5),
-    group = after_stat(x)
+    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 941be0cf..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))
@@ -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)
 )
diff --git a/man/conjugate.Rd b/man/conjugate.Rd
index 152fc8ab..897c7d96 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,15 +73,19 @@ 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{support}{Deprecated}
 }
 \value{
 A \link{conjugate-class} object with slots including:
@@ -173,6 +174,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
@@ -203,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
@@ -212,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
@@ -229,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
@@ -238,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
@@ -302,6 +309,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(
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 be127b45..85a13fa6 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, plot = TRUE,
-      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)
@@ -338,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.4252331, 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))