Skip to content

Commit

Permalink
issue 84
Browse files Browse the repository at this point in the history
  • Loading branch information
xidongdxi committed Jul 10, 2024
1 parent 97f0372 commit a07b150
Show file tree
Hide file tree
Showing 32 changed files with 385 additions and 224 deletions.
14 changes: 8 additions & 6 deletions DESCRIPTION
Original file line number Diff line number Diff line change
@@ -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", , "dong.xi1@gilead.com", role = c("aut", "cre")),
person("Ethan", "Brockmann", , "ethan.brockmann@atorusresearch.com", 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
Expand Down
5 changes: 5 additions & 0 deletions NAMESPACE
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand Down
5 changes: 5 additions & 0 deletions NEWS.md
Original file line number Diff line number Diff line change
Expand Up @@ -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
47 changes: 23 additions & 24 deletions R/adjust_p.R
Original file line number Diff line number Diff line change
Expand Up @@ -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`.
Expand All @@ -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
Expand All @@ -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
Expand All @@ -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)
Expand All @@ -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,
Expand Down Expand Up @@ -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)
Expand Down
125 changes: 26 additions & 99 deletions R/adjust_weights.R
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand All @@ -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
Expand All @@ -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)
Expand All @@ -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,
Expand Down Expand Up @@ -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)]
#'
#' graphicalMCP:::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)

Expand All @@ -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)
}
Loading

0 comments on commit a07b150

Please sign in to comment.