Conj dirichlet

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)
 #' )

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)) {

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])

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)

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

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)

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)

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)

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