version 0.1.0

DESCRIPTION

@@ -0,0 +1,21 @@
+Package: bmet
+Title: Bayesian Multigroup Equivalence Testing
+Version: 0.1.0
+Authors@R: 
+    person(given = "Tony",
+           family = "Pourmohamad",
+           role = c("aut", "cre"),
+           email = "tpourmohamad@gmail.com")
+Description: Calculates the necessary quantities to perform Bayesian multigroup equivalence testing. 
+             Currently the package includes the Bayesian models and equivalence criteria outlined in Pourmohamad and Lee (2023) 
+             <doi:10.1002/sta4.645>, but more models and equivalence testing features may be added over time.
+License: MIT + file LICENSE
+Encoding: UTF-8
+RoxygenNote: 7.2.3
+Imports: MASS, MCMCpack, stats, utils
+NeedsCompilation: no
+Packaged: 2023-12-21 17:20:10 UTC; pourmoht
+Author: Tony Pourmohamad [aut, cre]
+Maintainer: Tony Pourmohamad <tpourmohamad@gmail.com>
+Repository: CRAN
+Date/Publication: 2024-01-08 23:20:02 UTC

LICENSE

@@ -0,0 +1,2 @@
+YEAR: 2023
+COPYRIGHT HOLDER: bmet authors

MD5

@@ -0,0 +1,11 @@
+2c69eb69fd5a5e1565ca0edc009fff55 *DESCRIPTION
+8b14c967f5c34109b7142a3f27c9eec8 *LICENSE
+919c87472fc8877f9f15a7b28221cc73 *NAMESPACE
+a46a477fbae8b360fdd944273f0ea850 *NEWS.md
+04be9b1f523371b0509b3660ed09d3d1 *R/apc.R
+0a41c58bde0fb60d1e054e50bc5b870a *R/bet.R
+fc43c81034e532062455a65360de0c1a *R/betpp.R
+47d2f874a11f339b1c5ad4a8cf69c012 *man/apc.Rd
+ff47865d1c98b6dc552137c55051127a *man/bet.Rd
+6cca8098aa5a209923d740884be9357f *man/betpp.Rd
+1668db3e22a276e2e15465b387d31968 *man/figures/logo.png

NAMESPACE

@@ -0,0 +1,5 @@
+# Generated by roxygen2: do not edit by hand
+
+export(apc)
+export(bet)
+export(betpp)

NEWS.md

@@ -0,0 +1,3 @@
+# bmet 0.1.0
+
+* Initial CRAN submission.

R/apc.R

@@ -0,0 +1,46 @@
+#' All pairwise comparisons
+#'
+#' Function creates a contrast matrix for all pairwise comparisons
+#'
+#' @param ngroups A positive integer greater than 1 denoting the number of groups
+#' @param labs A vector of groups labels with length equal to \code{ngroups}. The default is set to \code{NULL}, and if used, the labels will be set to \code{1:length(ngroups)}.
+#' @return The function returns a matrix of all pairwise contrasts.
+#' @examples
+#' ### A contrast matrix based on all pairwise contrasts of 5 groups
+#' apc(5)
+#'
+#' ### Adding group labels
+#' apc(5, labs = paste("Group", 1:5, sep = " "))
+#'
+#' @export
+apc <- function(ngroups, labs = NULL) {
+  if(((floor(ngroups) - ngroups) != 0) | ngroups < 2){
+    stop("The number of groups must be a positive integer greater than 1.")
+  }else if(length(labs) != ngroups & !is.null(labs)){
+    stop("The number of groups labels must be the same as the number of groups.")
+  }else{
+    lev <- labs
+    lfm <- diag(ngroups)
+    nlev <- nrow(lfm)
+    rn <- rownames(lfm)
+    a <- attr(lfm, "grid")
+    if(is.null(lev)){
+      if(!is.null(a)){
+        lev <- apply(a, 1, paste, collapse = ":")
+      }else if(!is.null(rn)){
+        lev <- rn
+      }else{
+        lev <- as.character(1:nlev)
+      }
+    }
+    cbn <- utils::combn(seq_along(lev), 2)
+    M <- lfm[cbn[1, ], ] - lfm[cbn[2, ], ]
+    if(is.vector(M)){
+      dim(M) <- c(1, length(M))
+    }
+    rownames(M) <- paste(lev[cbn[1, ]], lev[cbn[2, ]], sep = "-")
+    return(M)
+  }
+}
+
+

R/bet.R

@@ -0,0 +1,83 @@
+#' Posterior based Bayesian multigroup equivalence testing
+#'
+#' Function provides the necessary tools to carry out Bayesian multigroup equivalence testing based on sampling of the posterior distribution.
+#' The function returns posterior samples of the average differences amongst groups, as well as posterior samples of group variances.
+#'
+#' @param values A vector of measurements sorted in the same order as the \code{groups} variable.
+#' @param groups A vector of groups labels corresponding to the individual measurements in the \code{groups} variable.
+#' @param em A c x 2 matrix of lower and upper equivalence margins. Here, c is the number of pairwise comparisons of interest.
+#' @param A A c x k matrix of pairwise contrasts. Here, k is the number of groups, i.e., \code{length(unique(groups))}.
+#' @param B A positive integer specifying the number of posterior samples to draw. By default \code{B} is set to 10000.
+#' @param test Setting this to anything other than "mean" tells the function to not calculate the posterior probability that the average
+#'             differences fall within the equivalence margins (applicable when testing equivalence based on something other than just
+#'             average differences).
+#' @return The function returns a list object containing the following:
+#' \itemize{
+#'   \item prob: The posterior probability that the average differences fall within the equivalence margins. Only returned if \code{test == "mean"}.
+#'   \item delta: A B x c matrix of posterior samples of the average difference for each pairwise comparison of interest.
+#'   \item sigma2: A B x k matrix of posterior samples of the variance for each group.
+#' }
+#' @references Pourmohamad, T. and Lee, H.K.H. (2023). Equivalence Testing for Multiple Groups. Stat, e645.
+#' @examples
+#' ### Multigroup equivalence test for A vs. B and A vs. C
+#' values <- rnorm(75)
+#' groups <- rep(LETTERS[1:3], each = 25)
+#'
+#' mad1 <- 0.65  # The equivalence margin for A vs. B
+#' mad2 <- 0.65  # The equivalence margin for A vs. C
+#' mads <- c(mad1, mad2)
+#' mads <- cbind(-mads, mads)
+#'
+#' A <- apc(3)
+#' A <- A[1:2, ]
+#'
+#' out <- bet(values, groups, mads, A, B = 10000)
+#'
+#' out$prob   # The posterior probability that the average
+#'            # differences fall within the equivalence margins
+#'
+#' @export
+bet <- function(values, groups, em, A, B = 10000, test = "mean"){
+
+  if(((floor(B) - B) != 0) | B < 1){
+    stop("The number of posterior samples, B, must be a positive integer.")
+  }else if(length(groups) != length(values)){
+    stop("The length of values must be the same as the length of groups.")
+  }else{
+    mads <- em
+    n.contrasts <- nrow(A)
+    group.labels <- unique(groups)
+    groups <- factor(groups, levels = group.labels)
+    m <- length(group.labels)
+    n <- tapply(values, groups, length)
+
+    E.mu <- tapply(values, groups, mean)
+
+    sigma2 <- matrix(NA, ncol = m, nrow = B)
+    delta <- matrix(NA, ncol = n.contrasts, nrow = B)
+
+    for(i in 1:B){
+      for(j in 1:m){
+        x <- values[groups == group.labels[j]]
+        sigma2[i, j] <- MCMCpack::rinvgamma(1, n[j] / 2 - 1 / 2, sum((x - E.mu[j])^2) / 2)
+      }
+
+      Cov.mu <- diag(sigma2[i,] / n)
+      delta[i,] <- MASS::mvrnorm(1, mu = A %*% E.mu, Sigma = A %*% Cov.mu %*% t(A))
+
+    }
+
+    colnames(delta) <- rownames(A)
+    colnames(sigma2) <- group.labels
+
+    if(test == "mean"){
+      prob <- sum(rowSums(t(t(delta) > mads[,1]) + t(t(delta) < mads[,2])) == 2 * n.contrasts) / B
+    }else{
+      prob <- NA
+    }
+
+    return(list(prob = prob, delta = delta, sigma2 = sigma2))
+
+  }
+
+}

R/betpp.R

@@ -0,0 +1,73 @@
+#' Posterior predictive based Bayesian multigroup equivalence testing
+#'
+#' Function provides the necessary tools to carry out Bayesian multigroup equivalence testing based on sampling of the posterior predictive distribution.
+#' The function returns posterior predictive samples of future differences amongst groups.
+#'
+#' @param values A vector of measurements sorted in the same order as the \code{groups} variable.
+#' @param groups A vector of groups labels corresponding to the individual measurements in the \code{groups} variable.
+#' @param em A c x 2 matrix of lower and upper equivalence margins. Here, c is the number of pairwise comparisons of interest.
+#' @param A A c x k matrix of pairwise contrasts. Here, k is the number of groups, i.e., \code{length(unique(groups))}.
+#' @param B A positive integer specifying the number of posterior predictive samples to draw. By default \code{B} is set to 10000.
+#' @return The function returns a list object containing the following:
+#' \itemize{
+#'   \item prob: The probability that future differences fall within the equivalence margins.
+#'   \item delta: A B x c matrix of posterior predictive samples of future differences for each pairwise comparison of interest.
+#' }
+#' @references Pourmohamad, T. and Lee, H.K.H. (2023). Equivalence Testing for Multiple Groups. Stat, e645.
+#' @examples
+#' ### Multigroup equivalence test for A vs. B and A vs. C
+#' values <- rnorm(75)
+#' groups <- rep(LETTERS[1:3], each = 25)
+#'
+#' mad1 <- 0.65  # The equivalence margin for A vs. B
+#' mad2 <- 0.65  # The equivalence margin for A vs. C
+#' mads <- c(mad1, mad2)
+#' mads <- cbind(-mads, mads)
+#'
+#' A <- apc(3)
+#' A <- A[1:2, ]
+#'
+#' out <- betpp(values, groups, mads, A, B = 10000)
+#'
+#' out$prob   # The probability that future differences
+#'            # fall within the equivalence margins
+#'
+#' @export
+betpp <- function(values, groups, em, A, B = 10000){
+
+  if(((floor(B) - B) != 0) | B < 1){
+    stop("The number of posterior samples, B, must be a positive integer.")
+  }else if(length(groups) != length(values)){
+    stop("The length of values must be the same as the length of groups.")
+  }else{
+    mads <- em
+    n.contrasts <- nrow(A)
+    group.labels <- unique(groups)
+    groups <- factor(groups, levels = group.labels)
+    m <- length(group.labels)
+    n <- tapply(values, groups, length)
+
+    E.mu <- tapply(values, groups, mean)
+
+    sigma2 <- matrix(NA, ncol = m, nrow = B)
+    mu <- matrix(NA, ncol = m, nrow = B)
+    pp <- matrix(NA, ncol = n.contrasts, nrow = B)
+
+    for(i in 1:B){
+      for(j in 1:m){
+        x <- values[groups == group.labels[j]]
+        sigma2[i, j] <- MCMCpack::rinvgamma(1, n[j] / 2 - 1 / 2, sum((x - E.mu[j])^2) / 2)
+        mu[i, j] <- stats::rnorm(1, E.mu[j], sigma2[i, j] / n[j])
+      }
+
+      Cov.mu <- diag(sigma2[i,])
+      pp[i,] <- MASS::mvrnorm(1, mu = A %*% mu[i,], Sigma = A %*% Cov.mu %*% t(A))
+
+    }
+
+    prob <- sum(rowSums(t(t(pp) >= mads[,1]) + t(t(pp) <= mads[,2])) == 2 * n.contrasts) / B
+    return(list(prob = prob, ppdraws = pp))
+
+  }
+
+}