Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Cran resubmission #85

Merged
merged 8 commits into from
Jul 11, 2024
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
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", , "[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
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)]
#'
#' 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
Loading