Add DBpower to CRANhaven, because archived on 2025-01-11 23:02:00 +0000

DBpower/DESCRIPTION

@@ -0,0 +1,25 @@
+Package: DBpower
+Type: Package
+Title: Finite Sample Power Calculations for Detection Boundary Tests
+Version: 0.1.0
+Authors@R: person("Ryan", "Sun", email = "[email protected]", role
+        = c("aut", "cre"))
+Description: Calculates lower bound on power, upper bound on power, and
+        exact power (small sets only) for detection boundary tests
+        (e.g. Berk-Jones, Generalized Berk-Jones, innovated Berk-Jones)
+        used in set-based inference studies. These detection boundary
+        tests are described in Sun et al., (2019)
+        <doi:10.1080/01621459.2019.1660170>.
+License: GPL-3
+Encoding: UTF-8
+RoxygenNote: 7.1.1
+Imports: dplyr, magrittr, stats, mvtnorm, combinat, GBJ, kit,
+Suggests: knitr, bindata, rmarkdown
+VignetteBuilder: knitr
+NeedsCompilation: no
+Packaged: 2022-02-09 17:30:51 UTC; rsun3
+Author: Ryan Sun [aut, cre]
+Maintainer: Ryan Sun <[email protected]>
+Repository: CRAN
+Date/Publication: 2022-02-10 19:10:05 UTC
+Additional_repositories: https://cranhaven.r-universe.dev

DBpower/MD5

@@ -0,0 +1,39 @@
+bc85c9e39794691832f878b02e554f39 *DESCRIPTION
+8e3848bc1ef1294d87682e8bf49a0348 *NAMESPACE
+7bcd6c688e47a6f4d4527c44003c8b1f *R/calc_b1.R
+b17984ce5225ef7588888c2f510276af *R/calc_b2.R
+794c7b692db809b2cd5303948efadb0e *R/calc_exact_power.R
+28179885d73fc344af92be5c99abff03 *R/set_BJ_bounds.R
+634228112a842ec7d1b39ead97e41661 *R/set_GBJ_bounds.R
+dc75a94e4e53bd135be1494473478690 *R/sim_b1.R
+5eb37d4fb962b4cf5074ec1369a3ac31 *R/sim_b2.R
+b7cb99381c8144f355b0ed508c72e1ce *R/sim_power_indiv.R
+f59356d83b403e8205ce9ccfeb1fe797 *R/sim_power_mvn.R
+f7ca976a55c8b7106165cadb12ec727f *build/vignette.rds
+b3ec09ab82f939d329f302146529ef24 *inst/doc/DBpower_tutorial.R
+cf566a65e6ff2f175c50fd64f780eedb *inst/doc/DBpower_tutorial.Rmd
+02db03103bd5af69d2d217ddf19c6f8e *inst/doc/DBpower_tutorial.html
+52b16f0a28dc28bac93d0733ef5801bd *man/BJ_func.Rd
+a7f2382ecc13d19d3a993112b04f626b *man/GBJ_bound1_root.Rd
+9ea0df5aa22b0d621a56fbe0e19da93d *man/calc_exact_power.Rd
+69d3e824bc1f41c46f61c92cc6f1ddd1 *man/calcb1.Rd
+56d3003b7da7bcb4c6bf5c13b986f240 *man/calcb2.Rd
+990fc80572291580926b23acb8014a1d *man/checkBoundsCross.Rd
+8aa15b16536180f535d41f3cdff4aea1 *man/checkBoundsLower1.Rd
+ad3ab94e32c1c5cb1ceaebba449b9662 *man/checkBoundsLower2.Rd
+77ad002d5e75b6d7c7c2766e9bd46571 *man/checkBoundsUpper1.Rd
+6b9dc65c22eda02932b69f44aed814cd *man/checkBoundsUpper2.Rd
+bb4ab610ee2f113df354cb92eb198b39 *man/createMj.Rd
+dec0d859119d511dec3155aae02ebee1 *man/createMjk.Rd
+e7a2082f427b36f7aa2c1104adc360bd *man/performIntegralLower1.Rd
+51894d20c330a27381d2bcf29f3d179c *man/performIntegralLower2.Rd
+d44ca8017b141fca73154e332599adc8 *man/performIntegralUpper1.Rd
+e55c359bdfb4129617fbe15607e5ceab *man/performIntegralUpper2.Rd
+55702e567129fc1d70f84edb6fe3ecfd *man/set_BJ_bounds.Rd
+7e54021accbd9c84d6256588a3829c0b *man/set_GBJ_bounds.Rd
+7c8d64ae2398903500ecebc28f641b1f *man/sim_b1.Rd
+11c70cf0e1fda1f266aa11f08f086220 *man/sim_b2.Rd
+feb42f37aef3ef22f777d6d13fb428a0 *man/sim_power_mvn.Rd
+c3ee4eb5570ea6d30c2ffc7779567191 *man/sim_stats_mef.Rd
+6ad5331cea18601240c902f19bd0d540 *man/sim_stats_mo.Rd
+cf566a65e6ff2f175c50fd64f780eedb *vignettes/DBpower_tutorial.Rmd

DBpower/NAMESPACE

@@ -0,0 +1,22 @@
+# Generated by roxygen2: do not edit by hand
+
+export(calc_exact_power)
+export(calcb1)
+export(calcb2)
+export(createMj)
+export(createMjk)
+export(performIntegralLower1)
+export(performIntegralLower2)
+export(performIntegralUpper1)
+export(performIntegralUpper2)
+export(set_BJ_bounds)
+export(set_GBJ_bounds)
+export(sim_b1)
+export(sim_b2)
+export(sim_power_mvn)
+export(sim_stats_mef)
+export(sim_stats_mo)
+importFrom(magrittr,"%>%")
+importFrom(stats,gaussian)
+importFrom(stats,glm)
+importFrom(stats,rnorm)

DBpower/R/calc_b1.R

@@ -0,0 +1,153 @@
+#' calc_b1.R
+#'
+#' Calculate lower bound or upper bound on power when considering only the largest
+#' test statistic in magnitude, i.e. only |Z|_(J) and not |Z|_(J-1).
+#'
+#' @param lower Boolean, whether to calculate lower bound.
+#' @param upper Boolean, whether to calculate upper bound.
+#' @param muVec Mean vector of test statistics under the alternative (assuming it's MVN).
+#' @param sigMat Covariance matrix of test statistics under the alternative (assuming it's MVN).
+#' @param bounds A J*1 vector of bounds on the magnitudes of the test statistics, where
+#' the first element is the bound for |Z|_(1) and the last element is the bound for |Z|_(J).
+#'
+#' @return A list with the elements:
+#' \item{allProbsLower}{J*1 vector of all components summed to calculate lower bound.}
+#' \item{lowerProb}{Lower bound.}
+#' \item{allProbsUpper}{J*1 vector of all components summed to calculate upper bound.}
+#' \item{upperProb}{Upper bound.}
+#'
+#' @export
+#' @examples
+#' myCov <- matrix(data=0.3, nrow=5, ncol=5)
+#' diag(myCov) <- 1
+#' myBounds <- set_GBJ_bounds(alpha = 0.01, J=5, sig_vec = myCov[lower.tri(myCov)])
+#' calcb1(muVec = c(1, 0, 0, 0, 0), sigMat = myCov, bounds=myBounds)
+#'
+calcb1 <- function(lower = TRUE, upper = FALSE, muVec, sigMat, bounds) {
+
+  # size of set
+  J <- length(muVec)
+
+  # lower and upper bounds on integration. See Sun, Shi, & Lin for expressions used.
+  lowerBounds1 <- rep(0, 2*J - 1)
+  upperBounds1 <- c(bounds[J], rep(Inf, 2*J - 2))
+  lowerBounds2 <- c(-bounds[J], rep(-Inf, 2*J - 2))
+  upperBounds2 <- rep(0, 2*J - 1)
+
+  # if calculating lower bound
+  if (lower) {
+    # calculate the probability/integral for each partition of the subspace we care about
+    allProbsLower <- sapply(X = 1:J, FUN = performIntegralLower1, muVec = muVec, sigMat = sigMat,
+                       lBounds1 = lowerBounds1, uBounds1 = upperBounds1, lBounds2 = lowerBounds2,
+                       uBounds2 = upperBounds2)
+    # only take the positive probabilities
+    lowerProb <- 1 - sum(allProbsLower[which(allProbsLower > 0)])
+  } else {
+    allProbsLower <- NA
+    lowerProb <- NA
+  }
+
+  # if calculating upper bound
+  if (upper) {
+    # add to the bounds the -b1 \leq x \leq b1 part
+    lowerExt <- rep(-bounds[1], J - 1)
+    upperExt <- rep(bounds[1], J - 1)
+    # calculate the probability/integral for each partition of the subspace we care about
+    allProbsUpper <- sapply(X = 1:J, FUN = performIntegralUpper1, muVec = muVec, sigMat = sigMat,
+                            lBounds1 = c(lowerBounds1, lowerExt), uBounds1 = c(upperBounds1, upperExt),
+                            lBounds2 = c(lowerBounds2, lowerExt), uBounds2 = c(upperBounds2, upperExt))
+    # only take the positive probabilities
+    upperProb <- 1 - sum(allProbsUpper[which(allProbsUpper > 0)])
+  } else {
+    allProbsUpper <- NA
+    upperProb <- NA
+  }
+
+  # return
+  return(list(allProbsLower = allProbsLower, lowerProb = lowerProb,
+              allProbsUpper = allProbsUpper, upperProb = upperProb))
+}
+
+
+#' Create the matrix that linearly transforms the vector of test statistics
+#' into a quantity amenable for pmvnorm.
+#'
+#' @param j The element of the vector that is the largest.
+#' @param size The length of the set.
+#'
+#' @return The transformation matrix of dimension (2J-1)*(2J-1)
+#'
+#' @export
+#' @examples
+#' createMj(j=3, size=5)
+createMj <- function(j, size) {
+  # the Zj - Zk part (referred to as Zm - Zj in paper)
+  topHalf <- diag(rep(-1, size))[-j, ]
+  topHalf[, j] <- 1
+  # the Zj + Zk part (referred to as Zm + Zj in paper)
+  bottomHalf <- abs(topHalf)
+
+  # put it together
+  Mj <- rbind(0, topHalf, bottomHalf)
+  # for bound on Zm
+  Mj[1, j] <- 1
+  return(Mj)
+}
+
+
+#' Apply this function over 1:J to calculate each portion of the integral
+#' we need for the lower bound.
+#'
+#' @param j Apply over this integer, the element that will be the largest in magnitude.
+#' @param muVec Mean vector of test statistics under the alternative (assuming it's MVN).
+#' @param sigMat Covariance matrix of test statistics under the alternative (assuming it's MVN).
+#' @param lBounds1 A 2J-1 vector of lower bounds for the first integral (see paper).
+#' @param uBounds1 A 2J-1 vector of upper bounds for the second integral (see paper).
+#' @param lBounds2 A 2J-1 vector of lower bounds for the first integral (see paper).
+#' @param uBounds2 A 2J-1 vector of upper bounds for the second integral (see paper).
+#'
+#' @return The value of the integration.
+#'
+#' @export
+performIntegralLower1 <- function(j, muVec, sigMat, lBounds1, uBounds1, lBounds2, uBounds2) {
+  # transform the mean and variance using the transformMat created above
+  transformMat <- createMj(j = j, size = length(muVec))
+  newMu <- as.numeric(transformMat %*% muVec)
+  newSig <- transformMat %*% sigMat %*% t(transformMat)
+  # calculate integral
+  intOut1 <- mvtnorm::pmvnorm(lower = lBounds1, upper = uBounds1, mean = newMu, sigma = newSig)
+  intOut2 <- mvtnorm::pmvnorm(lower = lBounds2, upper = uBounds2, mean = newMu, sigma = newSig)
+  return(intOut1[1] + intOut2[1])
+}
+
+
+
+#' Apply this function over 1:J to calculate each portion of the integral
+#' we need for the upper bound.
+#'
+#' @param j Apply over this integer, the element that will be the largest in magnitude.
+#' @param muVec Mean vector of test statistics under the alternative (assuming it's MVN).
+#' @param sigMat Covariance matrix of test statistics under the alternative (assuming it's MVN).
+#' @param lBounds1 A 3J-2 vector of lower bounds for the first integral (see paper), bounds will be longer than for performIntegralLower1.
+#' @param uBounds1 A 3J-2 vector of upper bounds for the second integral (see paper).
+#' @param lBounds2 A 3J-2 vector of lower bounds for the first integral (see paper).
+#' @param uBounds2 A 3J-2 vector of upper bounds for the second integral (see paper).
+#'
+#' @return The value of the integration.
+#'
+#' @export
+performIntegralUpper1 <- function(j, muVec, sigMat, lBounds1, uBounds1, lBounds2, uBounds2) {
+  # transform the mean and variance using the transformMat created above
+  transformMatTop <- createMj(j = j, size = length(muVec))
+  addBottom <- diag(rep(1, length(muVec)))[-j, ]
+  # for the upper bound, we need to add J-1 of the -b1 \leq Zj \leq b1 parts
+  transformMat <- rbind(transformMatTop, addBottom)
+
+  newMu <- as.numeric(transformMat %*% muVec)
+  newSig <- transformMat %*% sigMat  %*% t(transformMat)
+
+  # bounds will be larger than for performIntegralLower1
+  intOut1 <- mvtnorm::pmvnorm(lower = lBounds1, upper = uBounds1, mean = newMu, sigma = newSig)
+  intOut2 <- mvtnorm::pmvnorm(lower = lBounds2, upper = uBounds2, mean = newMu, sigma = newSig)
+  return(intOut1[1] + intOut2[1])
+}