Merge pull request #85 from Gilead-BioStats/CRAN-resubmission

Cran resubmission

DESCRIPTION

@@ -1,17 +1,19 @@
 Type: Package
 Package: graphicalMCP
 Title: Graphical Multiple Comparison Procedures
-Version: 0.2.3
+Version: 0.2.4
 Authors@R: c(
     person("Dong", "Xi", , "[email protected]", role = c("aut", "cre")),
     person("Ethan", "Brockmann", , "[email protected]", role = "aut"),
     person("Gilead Sciences, Inc.", role = c("cph", "fnd"))
   )
-Description: Graphical multiple comparison procedures (MCPs) control
-    the familywise error rate. This class includes many commonly used 
-    procedures as special cases. This package is a low-dependency
-    implementation of graphical MCPs which allow mixed types of tests.
-    It also includes power simulations and visualization of graphical MCPs.
+Description: Multiple comparison procedures (MCPs) control the familywise error
+    rate in clinical trials. Graphical MCPs include many commonly used
+    procedures as special cases; see Bretz et al. (2011)
+    <doi:10.1002/bimj.201000239>, Lu (2016) <doi:10.1002/sim.6985>, and Xi et
+    al. (2017) <doi:10.1002/bimj.201600233>. This package is a low-dependency
+    implementation of graphical MCPs which allow mixed types of tests. It also
+    includes power simulations and visualization of graphical MCPs.
 License: Apache License (>= 2)
 URL: https://gilead-biostats.github.io/graphicalMCP/
 BugReports: https://github.com/Gilead-BioStats/graphicalMCP/issues

NAMESPACE

@@ -10,6 +10,11 @@ S3method(print,graph_report)
 S3method(print,initial_graph)
 S3method(print,power_report)
 S3method(print,updated_graph)
+export(adjust_p_bonferroni)
+export(adjust_p_parametric)
+export(adjust_p_simes)
+export(adjust_weights_parametric)
+export(adjust_weights_simes)
 export(as_graphMCP)
 export(as_igraph)
 export(as_initial_graph)

NEWS.md

@@ -31,3 +31,8 @@
 
 * Included cran-comments.ms in .Rbuildignore
 * Resubmission for first CRAN release
+
+# graphicalMCP 0.2.4
+
+* Updated documentation according to issue #84
+* Resubmission for first CRAN release

R/adjust_p.R

@@ -6,10 +6,10 @@
 #' The intersection hypothesis can be rejected if its adjusted p-value is less
 #' than or equal to \eqn{\alpha}. Currently, there are three test types
 #' supported:
-#' * Bonferroni tests for `graphicalMCP:::adjust_p_bonferroni()`,
-#' * Parametric tests for `graphicalMCP:::adjust_p_parametric()`,
+#' * Bonferroni tests for [adjust_p_bonferroni()],
+#' * Parametric tests for [adjust_p_parametric()],
 #'     - Note that one-sided tests are required for parametric tests.
-#' * Simes tests for `graphicalMCP:::adjust_p_simes()`.
+#' * Simes tests for [adjust_p_simes()].
 #'
 #' @param p A numeric vector of p-values (unadjusted, raw), whose values should
 #'   be between 0 & 1. The length should match the length of `hypotheses`.
@@ -18,7 +18,7 @@
 #'   `p`. The sum of hypothesis weights should not exceed 1.
 #' @param test_corr (Optional) A numeric matrix of correlations between test
 #'   statistics, which is needed to perform parametric tests using
-#'   `graphicalMCP:::adjust_p_parametric()`. The number of rows and columns of
+#'   [adjust_p_parametric()]. The number of rows and columns of
 #'   this correlation matrix should match the length of `p`.
 #' @param maxpts (Optional) An integer scalar for the maximum number of function
 #'   values, which is needed to perform parametric tests using the
@@ -33,13 +33,13 @@
 #' @return A single adjusted p-value for the intersection hypothesis.
 #'
 #' @seealso
-#'   `graphicalMCP:::adjust_weights_parametric()` for adjusted hypothesis
-#'   weights using parametric tests, `graphicalMCP:::adjust_weights_simes()`
-#'   for adjusted hypothesis weights using Simes tests.
+#'   [adjust_weights_parametric()] for adjusted hypothesis weights using
+#'   parametric tests, [adjust_weights_simes()] for adjusted hypothesis weights
+#'   using Simes tests.
 #'
 #' @rdname adjust_p
 #'
-#' @keywords internal
+#' @export
 #'
 #' @references
 #'   Bretz, F., Maurer, W., Brannath, W., and Posch, M. (2009). A graphical
@@ -54,22 +54,9 @@
 #'   \emph{Biometrical Journal}, 59(5), 918-931.
 #'
 #' @examples
-#' set.seed(1234)
-#'
 #' hypotheses <- c(H1 = 0.5, H2 = 0.25, H3 = 0.25)
 #' p <- c(0.019, 0.025, 0.05)
-#'
-#' # Bonferroni test
-#' graphicalMCP:::adjust_p_bonferroni(p, hypotheses)
-#'
-#' # Simes test
-#' graphicalMCP:::adjust_p_simes(p, hypotheses)
-#'
-#' # Parametric test
-#' # Using the `mvtnorm::GenzBretz` algorithm
-#' corr <- matrix(0.5, nrow = 3, ncol = 3)
-#' diag(corr) <- 1
-#' graphicalMCP:::adjust_p_parametric(p, hypotheses, corr)
+#' adjust_p_bonferroni(p, hypotheses)
 adjust_p_bonferroni <- function(p, hypotheses) {
   if (sum(hypotheses) == 0) {
     return(Inf)
@@ -84,7 +71,15 @@ adjust_p_bonferroni <- function(p, hypotheses) {
 }
 
 #' @rdname adjust_p
-#' @keywords internal
+#' @export
+#' @examples
+#' set.seed(1234)
+#' hypotheses <- c(H1 = 0.5, H2 = 0.25, H3 = 0.25)
+#' p <- c(0.019, 0.025, 0.05)
+#' # Using the `mvtnorm::GenzBretz` algorithm
+#' corr <- matrix(0.5, nrow = 3, ncol = 3)
+#' diag(corr) <- 1
+#' adjust_p_parametric(p, hypotheses, corr)
 adjust_p_parametric <- function(p,
                                 hypotheses,
                                 test_corr = NULL,
@@ -121,7 +116,11 @@ adjust_p_parametric <- function(p,
 }
 
 #' @rdname adjust_p
-#' @keywords internal
+#' @export
+#' @examples
+#' hypotheses <- c(H1 = 0.5, H2 = 0.25, H3 = 0.25)
+#' p <- c(0.019, 0.025, 0.05)
+#' adjust_p_simes(p, hypotheses)
 adjust_p_simes <- function(p, hypotheses) {
   if (sum(hypotheses) == 0) {
     return(Inf)

R/adjust_weights.R

@@ -6,31 +6,27 @@
 #' hypothesis weights times \eqn{\alpha}. For Bonferroni tests, their adjusted
 #' hypothesis weights are their hypothesis weights of the intersection
 #' hypothesis. Additional adjustment is needed for parametric and Simes tests:
-#' * Parametric tests for `graphicalMCP:::adjust_weights_parametric()`,
+#' * Parametric tests for [adjust_weights_parametric()],
 #'     - Note that one-sided tests are required for parametric tests.
-#' * Simes tests for `graphicalMCP:::adjust_weights_simes()`.
+#' * Simes tests for [adjust_weights_simes()].
 #'
 #' @param matrix_weights (Optional) A matrix of hypothesis weights of all
 #'   intersection hypotheses. This can be obtained as the second half of columns
 #'   from the output of [graph_generate_weights()].
 #' @param matrix_intersections (Optional) A matrix of hypothesis indicators of
 #'   all intersection hypotheses. This can be obtained as the first half of
 #'   columns from the output of [graph_generate_weights()].
-#' @param x The root to solve for with `stats::uniroot()`.
 #' @param alpha (Optional) A numeric value of the overall significance level,
 #'   which should be between 0 & 1. The default is 0.025 for one-sided
 #'   hypothesis testing problems; another common choice is 0.05 for two-sided
 #'   hypothesis testing problems. Note when parametric tests are used, only
 #'   one-sided tests are supported.
 #' @param p (Optional) A numeric vector of p-values (unadjusted, raw), whose
-#'   values should be between 0 & 1. The length should match the length of
-#'   `hypotheses`.
-#' @param hypotheses (Optional) A numeric vector of hypothesis weights. Must be
-#'   a vector of values between 0 & 1 (inclusive). The length should match the
-#'   length of `p`. The sum of hypothesis weights should not exceed 1.
+#'   values should be between 0 & 1. The length should match the number of
+#'   columns of `matrix_weights`.
 #' @param test_corr (Optional) A numeric matrix of correlations between test
 #'   statistics, which is needed to perform parametric tests using
-#'   `graphicalMCP:::adjust_p_parametric()`. The number of rows and columns of
+#'   [adjust_weights_parametric()]. The number of rows and columns of
 #'   this correlation matrix should match the length of `p`.
 #' @param test_groups (Optional) A list of numeric vectors specifying hypotheses
 #'   to test together. Grouping is needed to correctly perform Simes and
@@ -46,28 +42,20 @@
 #'   `mvtnorm::GenzBretz` algorithm. The default is 0.
 #'
 #' @return
-#' * `graphicalMCP:::adjust_weights_parametric()` returns a matrix with the same
+#' * [adjust_weights_parametric()] returns a matrix with the same
 #'   dimensions as `matrix_weights`, whose hypothesis weights have been
 #'   adjusted according to parametric tests.
-#' * `graphicalMCP:::adjust_weights_simes()` returns a matrix with the same
+#' * [adjust_weights_simes()] returns a matrix with the same
 #'   dimensions as `matrix_weights`, whose hypothesis weights have been
 #'   adjusted according to Simes tests.
-#' * `graphicalMCP:::c_value_function()` returns the difference between
-#'   \eqn{\alpha} and the Type I error of the parametric test with the \eqn{c}
-#'   value of `x`, adjusted for the correlation between test statistics using
-#'   parametric tests based on equation (6) of Xi et al. (2017).
-#' * `graphicalMCP:::solve_c_parametric()` returns the c value adjusted for the
-#'   correlation between test statistics using parametric tests based on
-#'   equation (6) of Xi et al. (2017).
 #'
 #' @seealso
-#'   `graphicalMCP:::adjust_p_parametric()` for adjusted p-values using
-#'   parametric tests, `graphicalMCP:::adjust_p_simes()` for adjusted p-values
-#'   using Simes tests.
+#'   [adjust_p_parametric()] for adjusted p-values using parametric tests,
+#'   [adjust_p_simes()] for adjusted p-values using Simes tests.
 #'
 #' @rdname adjust_weights
 #'
-#' @keywords internal
+#' @export
 #'
 #' @references
 #'   Lu, K. (2016). Graphical approaches using a Bonferroni mixture of weighted
@@ -78,7 +66,6 @@
 #'   \emph{Biometrical Journal}, 59(5), 918-931.
 #'
 #' @examples
-#' set.seed(1234)
 #' alpha <- 0.025
 #' p <- c(0.018, 0.01, 0.105, 0.006)
 #' num_hyps <- length(p)
@@ -88,25 +75,13 @@
 #' matrix_weights <- weighting_strategy[, -seq_len(num_hyps)]
 #'
 #' set.seed(1234)
-#' graphicalMCP:::adjust_weights_parametric(
+#' adjust_weights_parametric(
 #'   matrix_weights = matrix_weights,
 #'   matrix_intersections = matrix_intersections,
 #'   test_corr = diag(4),
 #'   alpha = alpha,
 #'   test_groups = list(1:4)
 #' )
-#'
-#' graphicalMCP:::adjust_weights_simes(
-#'   matrix_weights = matrix_weights,
-#'   p = p,
-#'   test_groups = list(1:4)
-#' )
-#'
-#' graphicalMCP:::solve_c_parametric(
-#'   hypotheses = matrix_weights[1, ],
-#'   test_corr = diag(4),
-#'   alpha = alpha
-#' )
 adjust_weights_parametric <- function(matrix_weights,
                                       matrix_intersections,
                                       test_corr,
@@ -146,7 +121,21 @@ adjust_weights_parametric <- function(matrix_weights,
 }
 
 #' @rdname adjust_weights
-#' @keywords internal
+#' @export
+#' @examples
+#' alpha <- 0.025
+#' p <- c(0.018, 0.01, 0.105, 0.006)
+#' num_hyps <- length(p)
+#' g <- bonferroni_holm(rep(1 / 4, 4))
+#' weighting_strategy <- graph_generate_weights(g)
+#' matrix_intersections <- weighting_strategy[, seq_len(num_hyps)]
+#' matrix_weights <- weighting_strategy[, -seq_len(num_hyps)]
+#'
+#' adjust_weights_simes(
+#'   matrix_weights = matrix_weights,
+#'   p = p,
+#'   test_groups = list(1:4)
+#' )
 adjust_weights_simes <- function(matrix_weights, p, test_groups) {
   ordered_p <- order(p)
 
@@ -162,65 +151,3 @@ adjust_weights_simes <- function(matrix_weights, p, test_groups) {
 
   do.call(cbind, group_adjusted_weights)
 }
-
-#' @rdname adjust_weights
-#' @keywords internal
-c_value_function <- function(x,
-                             hypotheses,
-                             test_corr,
-                             alpha,
-                             maxpts = 25000,
-                             abseps = 1e-6,
-                             releps = 0) {
-  hyps_nonzero <- which(hypotheses > 0)
-  z <- stats::qnorm(x * hypotheses[hyps_nonzero] * alpha, lower.tail = FALSE)
-  y <- ifelse(
-    length(z) == 1,
-    stats::pnorm(z, lower.tail = FALSE)[[1]],
-    1 - mvtnorm::pmvnorm(
-      lower = -Inf,
-      upper = z,
-      corr = test_corr[hyps_nonzero, hyps_nonzero, drop = FALSE],
-      algorithm = mvtnorm::GenzBretz(
-        maxpts = maxpts,
-        abseps = abseps,
-        releps = releps
-      )
-    )[[1]]
-  )
-
-  y - alpha * sum(hypotheses)
-}
-
-#' @rdname adjust_weights
-#' @keywords internal
-solve_c_parametric <- function(hypotheses,
-                               test_corr,
-                               alpha,
-                               maxpts = 25000,
-                               abseps = 1e-6,
-                               releps = 0) {
-  num_hyps <- seq_along(hypotheses)
-  c_value <- ifelse(
-    length(num_hyps) == 1 || sum(hypotheses) == 0,
-    1,
-    stats::uniroot(
-      c_value_function,
-      lower = 0, # Why is this not -Inf? Ohhhh because c_value >= 1
-      # upper > 40 errors when w[i] ~= 1 && w[j] = epsilon
-      # upper = 2 errors when w = c(.5, .5) && all(test_corr == 1)
-      # furthermore, even under perfect correlation & with balanced weights, the
-      # c_function_to_solve does not exceed `length(hypotheses)`
-      upper = length(hypotheses) + 1,
-      hypotheses = hypotheses,
-      test_corr = test_corr,
-      alpha = alpha,
-      maxpts = maxpts,
-      abseps = abseps,
-      releps = releps
-    )$root
-  )
-
-  # Occasionally has floating point differences
-  round(c_value, 10)
-}