diff --git a/DESCRIPTION b/DESCRIPTION index c5433bc..72fb736 100644 --- a/DESCRIPTION +++ b/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", , "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) + , Lu (2016) , and Xi et + al. (2017) . 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 diff --git a/NAMESPACE b/NAMESPACE index 0420fbf..ecf1472 100644 --- a/NAMESPACE +++ b/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) diff --git a/NEWS.md b/NEWS.md index d2582cb..c724408 100644 --- a/NEWS.md +++ b/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 diff --git a/R/adjust_p.R b/R/adjust_p.R index c3f0749..8077144 100644 --- a/R/adjust_p.R +++ b/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) diff --git a/R/adjust_weights.R b/R/adjust_weights.R index ca64bfe..a12f3df 100644 --- a/R/adjust_weights.R +++ b/R/adjust_weights.R @@ -6,9 +6,9 @@ #' 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 @@ -16,21 +16,17 @@ #' @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)] +#' +#' 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) @@ -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) -} diff --git a/R/adjust_weights_parametric_util.R b/R/adjust_weights_parametric_util.R new file mode 100644 index 0000000..a67c3f3 --- /dev/null +++ b/R/adjust_weights_parametric_util.R @@ -0,0 +1,120 @@ +#' Calculate adjusted hypothesis weights for parametric tests +#' +#' @description +#' An intersection hypothesis can be rejected if its p-values are less than or +#' equal to their adjusted significance levels, which are their adjusted +#' 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 tests: +#' * Parametric tests for [adjust_weights_parametric()], +#' - Note that one-sided tests are required for parametric tests. +#' +#' @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 hypotheses 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. +#' @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 test_corr (Optional) A numeric matrix of correlations between test +#' statistics, which is needed to perform parametric tests using +#' [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 +#' parametric tests. +#' @param maxpts (Optional) An integer scalar for the maximum number of function +#' values, which is needed to perform parametric tests using the +#' `mvtnorm::GenzBretz` algorithm. The default is 25000. +#' @param abseps (Optional) A numeric scalar for the absolute error tolerance, +#' which is needed to perform parametric tests using the `mvtnorm::GenzBretz` +#' algorithm. The default is 1e-6. +#' @param releps (Optional) A numeric scalar for the relative error tolerance +#' as double, which is needed to perform parametric tests using the +#' `mvtnorm::GenzBretz` algorithm. The default is 0. +#' +#' @return +#' * `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). +#' * `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 +#' [adjust_weights_parametric()] for adjusted hypothesis weights using +#' parametric tests. +#' +#' @rdname adjust_weights_parametric_util +#' +#' @keywords internal +#' +#' @references +#' Xi, D., Glimm, E., Maurer, W., and Bretz, F. (2017). A unified framework +#' for weighted parametric multiple test procedures. +#' \emph{Biometrical Journal}, 59(5), 918-931. +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_parametric_util +#' @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) +} diff --git a/R/edge_pairs.R b/R/edge_pairs.R index 3b81acd..c871ee6 100644 --- a/R/edge_pairs.R +++ b/R/edge_pairs.R @@ -13,8 +13,6 @@ #' #' @keywords internal #' -#' @examples -#' graphicalMCP:::edge_pairs(bonferroni_holm(hypotheses = rep(1 / 3, 3))) edge_pairs <- function(graph) { g_names <- names(graph$hypotheses) diff --git a/R/graph_calculate_power.R b/R/graph_calculate_power.R index 92fb8fb..df8d5d9 100644 --- a/R/graph_calculate_power.R +++ b/R/graph_calculate_power.R @@ -132,11 +132,12 @@ #' ) #' set.seed(1234) #' # Bonferroni tests +#' # Reduce the number of simulations to save time for package compilation #' power_output <- graph_calculate_power( #' g, #' alpha, #' sim_corr = corr, -#' sim_n = 1e5, +#' sim_n = 1e2, #' power_marginal = marginal_power, #' sim_success = success_fns #' ) @@ -150,7 +151,7 @@ #' test_groups = list(1:2, 3:4), #' test_types = c("parametric", "simes"), #' test_corr = list(corr1, NA), -#' sim_n = 1e3, +#' sim_n = 1e2, #' sim_success = list( #' function(.) .[1] || .[2], #' function(.) .[1] && .[2] diff --git a/R/graph_create.R b/R/graph_create.R index 2a48e00..33e2b00 100644 --- a/R/graph_create.R +++ b/R/graph_create.R @@ -101,7 +101,7 @@ #' H3 = c(0, 1, 0, 0), #' H4 = c(1, 0, 0, 0) #' ) -#' # g <- graph_create(hypotheses, transitions) +#' \dontrun{g <- graph_create(hypotheses, transitions)} #' #' # When names are not specified, hypotheses are numbered sequentially as #' # H1, H2, ... diff --git a/R/plot.initial_graph.R b/R/plot.initial_graph.R index 65d2540..061fbd6 100644 --- a/R/plot.initial_graph.R +++ b/R/plot.initial_graph.R @@ -38,7 +38,8 @@ #' Defaults to all 0, since igraph plots tend to have large margins. It is #' passed directly to [graphics::par()] (`mar`). #' -#' @return NULL, after plotting the initial graph. +#' @return An object x of class `initial_graph`, after plotting the initial +#' graph. #' #' @section Customization of graphs: #' There are a few values for [igraph::plot.igraph()] that get their defaults @@ -232,4 +233,6 @@ plot.initial_graph <- function(x, edge.arrow.width = 1, asp = 0 ) + + invisible(x) } diff --git a/R/plot.updated_graph.R b/R/plot.updated_graph.R index 0dc7885..7041749 100644 --- a/R/plot.updated_graph.R +++ b/R/plot.updated_graph.R @@ -8,7 +8,8 @@ #' @param x An object of class `updated_graph` to plot. #' @inheritDotParams plot.initial_graph #' -#' @return NULL, after plotting the updated graph. +#' @return An object x of class `updated_graph`, after plotting the updated +#' graph. #' #' @seealso #' [plot.initial_graph()] for the plot method for the initial graph. @@ -49,4 +50,6 @@ plot.updated_graph <- function(x, ...) { v_colors[x$deleted] <- "#cccccc" plot(x$updated_graph, vertex.color = v_colors, ...) + + invisible(x) } diff --git a/R/print.graph_report.R b/R/print.graph_report.R index f0070b8..9bd9b30 100644 --- a/R/print.graph_report.R +++ b/R/print.graph_report.R @@ -12,6 +12,9 @@ #' @param rows An integer scalar indicating how many rows of detailed test #' results to print. #' +#' @return An object x of class `graph_report`, after printing the report of +#' conducting a graphical multiple comparison procedure. +#' #' @rdname print.graph_report #' #' @export diff --git a/R/print.initial_graph.R b/R/print.initial_graph.R index fdbb045..2b32367 100644 --- a/R/print.initial_graph.R +++ b/R/print.initial_graph.R @@ -10,7 +10,8 @@ #' to to display. #' @param indent An integer scalar indicating how many spaces to indent results. #' -#' @return NULL, after printing the initial graph. +#' @return An object x of class `initial_graph`, after printing the initial +#' graph. #' #' @seealso #' [print.updated_graph()] for the print method for the updated graph after @@ -92,4 +93,6 @@ print.initial_graph <- function(x, transitions_text <- data.frame(df_trn, check.names = FALSE) print(transitions_text, row.names = FALSE) + + invisible(x) } diff --git a/R/print.power_report.R b/R/print.power_report.R index a9bfc80..cc1028e 100644 --- a/R/print.power_report.R +++ b/R/print.power_report.R @@ -7,6 +7,10 @@ #' @param x An object of the class `power_report` to print #' @inheritParams print.graph_report #' +#' @return An object x of the class `power_report`, after printing the report of +#' conducting power simulations based on a graphical multiple comparison +#' procedure. +#' #' @rdname print.power_report #' #' @export diff --git a/R/print.updated_graph.R b/R/print.updated_graph.R index 6ab737b..3b1f809 100644 --- a/R/print.updated_graph.R +++ b/R/print.updated_graph.R @@ -11,7 +11,8 @@ #' to to display. #' @param indent An integer scalar indicating how many spaces to indent results. #' -#' @return NULL, after printing the updated graph. +#' @return An object x of the class `updated_graph`, after printing the updated +#' graph. #' #' @seealso #' [print.initial_graph()] for the print method for the initial graph. @@ -107,4 +108,6 @@ print.updated_graph <- function(x, ..., precision = 6, indent = 2) { ) cat("\n") } + + invisible(x) } diff --git a/_pkgdown.yml b/_pkgdown.yml index cd7fbb8..b9ea367 100644 --- a/_pkgdown.yml +++ b/_pkgdown.yml @@ -30,6 +30,11 @@ reference: - graph_test_shortcut - print.graph_report - graph_rejection_orderings + - adjust_p_bonferroni + - adjust_p_parametric + - adjust_p_simes + - adjust_weights_parametric + - adjust_weights_simes - title: Power simulation - desc: Functions for performing power simulations - contents: diff --git a/codemeta.json b/codemeta.json index 486c0b7..e7b25b9 100644 --- a/codemeta.json +++ b/codemeta.json @@ -2,12 +2,12 @@ "@context": "https://doi.org/10.5063/schema/codemeta-2.0", "@type": "SoftwareSourceCode", "identifier": "graphicalMCP", - "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) , Lu (2016) , and Xi et al. (2017) . 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.", "name": "graphicalMCP: Graphical Multiple Comparison Procedures", "codeRepository": "https://github.com/Gilead-BioStats/graphicalMCP", "issueTracker": "https://github.com/Gilead-BioStats/graphicalMCP/issues", "license": "Apache License 2", - "version": "0.2.3", + "version": "0.2.4", "programmingLanguage": { "@type": "ComputerLanguage", "name": "R", @@ -312,7 +312,7 @@ }, "SystemRequirements": null }, - "fileSize": "713.316KB", + "fileSize": "727.413KB", "citation": [ { "@type": "SoftwareSourceCode", diff --git a/cran-comments.md b/cran-comments.md index 7b2f946..f76a6cb 100644 --- a/cran-comments.md +++ b/cran-comments.md @@ -1,3 +1,14 @@ +# Version 0.2.4 + +- This is the fourth submission to CRAN. + + Updated documentation according to issue #84 + + There is one note regarding the spelling of "MCPs" and "familywise". This is + to confirm that they are correct. + +## R CMD check results + +0 errors | 0 warnings | 0 notes + # Version 0.2.3 - This is the third submission to CRAN. diff --git a/inst/CITATION b/inst/CITATION index bffa8be..8c720ca 100644 --- a/inst/CITATION +++ b/inst/CITATION @@ -4,7 +4,7 @@ bibentry( "manual", title = "{graphicalMCP}: Graphical multiple comparison procedures", author = as.person("Dong Xi and Ethan Brockmann"), - edition = "version 0.2.3", + edition = "version 0.2.4", year = 2024, url = "https://github.com/Gilead-BioStats/graphicalMCP" ) diff --git a/man/adjust_p.Rd b/man/adjust_p.Rd index 08a5a5f..6cbcac4 100644 --- a/man/adjust_p.Rd +++ b/man/adjust_p.Rd @@ -29,7 +29,7 @@ values between 0 & 1 (inclusive). The length should match the length of \item{test_corr}{(Optional) A numeric matrix of correlations between test statistics, which is needed to perform parametric tests using -\code{graphicalMCP:::adjust_p_parametric()}. The number of rows and columns of +\code{\link[=adjust_p_parametric]{adjust_p_parametric()}}. The number of rows and columns of this correlation matrix should match the length of \code{p}.} \item{maxpts}{(Optional) An integer scalar for the maximum number of function @@ -54,31 +54,28 @@ 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: \itemize{ -\item Bonferroni tests for \code{graphicalMCP:::adjust_p_bonferroni()}, -\item Parametric tests for \code{graphicalMCP:::adjust_p_parametric()}, +\item Bonferroni tests for \code{\link[=adjust_p_bonferroni]{adjust_p_bonferroni()}}, +\item Parametric tests for \code{\link[=adjust_p_parametric]{adjust_p_parametric()}}, \itemize{ \item Note that one-sided tests are required for parametric tests. } -\item Simes tests for \code{graphicalMCP:::adjust_p_simes()}. +\item Simes tests for \code{\link[=adjust_p_simes]{adjust_p_simes()}}. } } \examples{ +hypotheses <- c(H1 = 0.5, H2 = 0.25, H3 = 0.25) +p <- c(0.019, 0.025, 0.05) +adjust_p_bonferroni(p, hypotheses) 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_parametric(p, hypotheses, corr) +hypotheses <- c(H1 = 0.5, H2 = 0.25, H3 = 0.25) +p <- c(0.019, 0.025, 0.05) +adjust_p_simes(p, hypotheses) } \references{ Bretz, F., Maurer, W., Brannath, W., and Posch, M. (2009). A graphical @@ -93,8 +90,7 @@ for weighted parametric multiple test procedures. \emph{Biometrical Journal}, 59(5), 918-931. } \seealso{ -\code{graphicalMCP:::adjust_weights_parametric()} for adjusted hypothesis -weights using parametric tests, \code{graphicalMCP:::adjust_weights_simes()} -for adjusted hypothesis weights using Simes tests. +\code{\link[=adjust_weights_parametric]{adjust_weights_parametric()}} for adjusted hypothesis weights using +parametric tests, \code{\link[=adjust_weights_simes]{adjust_weights_simes()}} for adjusted hypothesis weights +using Simes tests. } -\keyword{internal} diff --git a/man/adjust_weights.Rd b/man/adjust_weights.Rd index d5461aa..ca3ffdd 100644 --- a/man/adjust_weights.Rd +++ b/man/adjust_weights.Rd @@ -3,8 +3,6 @@ \name{adjust_weights_parametric} \alias{adjust_weights_parametric} \alias{adjust_weights_simes} -\alias{c_value_function} -\alias{solve_c_parametric} \title{Calculate adjusted hypothesis weights} \usage{ adjust_weights_parametric( @@ -19,25 +17,6 @@ adjust_weights_parametric( ) adjust_weights_simes(matrix_weights, p, test_groups) - -c_value_function( - x, - hypotheses, - test_corr, - alpha, - maxpts = 25000, - abseps = 1e-06, - releps = 0 -) - -solve_c_parametric( - hypotheses, - test_corr, - alpha, - maxpts = 25000, - abseps = 1e-06, - releps = 0 -) } \arguments{ \item{matrix_weights}{(Optional) A matrix of hypothesis weights of all @@ -50,7 +29,7 @@ columns from the output of \code{\link[=graph_generate_weights]{graph_generate_w \item{test_corr}{(Optional) A numeric matrix of correlations between test statistics, which is needed to perform parametric tests using -\code{graphicalMCP:::adjust_p_parametric()}. The number of rows and columns of +\code{\link[=adjust_weights_parametric]{adjust_weights_parametric()}}. The number of rows and columns of this correlation matrix should match the length of \code{p}.} \item{alpha}{(Optional) A numeric value of the overall significance level, @@ -76,30 +55,17 @@ as double, which is needed to perform parametric tests using the \code{mvtnorm::GenzBretz} algorithm. The default is 0.} \item{p}{(Optional) A numeric vector of p-values (unadjusted, raw), whose -values should be between 0 & 1. The length should match the length of -\code{hypotheses}.} - -\item{x}{The root to solve for with \code{stats::uniroot()}.} - -\item{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 \code{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 \code{matrix_weights}.} } \value{ \itemize{ -\item \code{graphicalMCP:::adjust_weights_parametric()} returns a matrix with the same +\item \code{\link[=adjust_weights_parametric]{adjust_weights_parametric()}} returns a matrix with the same dimensions as \code{matrix_weights}, whose hypothesis weights have been adjusted according to parametric tests. -\item \code{graphicalMCP:::adjust_weights_simes()} returns a matrix with the same +\item \code{\link[=adjust_weights_simes]{adjust_weights_simes()}} returns a matrix with the same dimensions as \code{matrix_weights}, whose hypothesis weights have been adjusted according to Simes tests. -\item \code{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 \code{x}, adjusted for the correlation between test statistics using -parametric tests based on equation (6) of Xi et al. (2017). -\item \code{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). } } \description{ @@ -109,15 +75,14 @@ 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: \itemize{ -\item Parametric tests for \code{graphicalMCP:::adjust_weights_parametric()}, +\item Parametric tests for \code{\link[=adjust_weights_parametric]{adjust_weights_parametric()}}, \itemize{ \item Note that one-sided tests are required for parametric tests. } -\item Simes tests for \code{graphicalMCP:::adjust_weights_simes()}. +\item Simes tests for \code{\link[=adjust_weights_simes]{adjust_weights_simes()}}. } } \examples{ -set.seed(1234) alpha <- 0.025 p <- c(0.018, 0.01, 0.105, 0.006) num_hyps <- length(p) @@ -127,25 +92,26 @@ matrix_intersections <- weighting_strategy[, seq_len(num_hyps)] 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) ) +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) ) - -graphicalMCP:::solve_c_parametric( - hypotheses = matrix_weights[1, ], - test_corr = diag(4), - alpha = alpha -) } \references{ Lu, K. (2016). Graphical approaches using a Bonferroni mixture of weighted @@ -156,8 +122,6 @@ for weighted parametric multiple test procedures. \emph{Biometrical Journal}, 59(5), 918-931. } \seealso{ -\code{graphicalMCP:::adjust_p_parametric()} for adjusted p-values using -parametric tests, \code{graphicalMCP:::adjust_p_simes()} for adjusted p-values -using Simes tests. +\code{\link[=adjust_p_parametric]{adjust_p_parametric()}} for adjusted p-values using parametric tests, +\code{\link[=adjust_p_simes]{adjust_p_simes()}} for adjusted p-values using Simes tests. } -\keyword{internal} diff --git a/man/adjust_weights_parametric_util.Rd b/man/adjust_weights_parametric_util.Rd new file mode 100644 index 0000000..b4ef3c9 --- /dev/null +++ b/man/adjust_weights_parametric_util.Rd @@ -0,0 +1,98 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/adjust_weights_parametric_util.R +\name{c_value_function} +\alias{c_value_function} +\alias{solve_c_parametric} +\title{Calculate adjusted hypothesis weights for parametric tests} +\usage{ +c_value_function( + x, + hypotheses, + test_corr, + alpha, + maxpts = 25000, + abseps = 1e-06, + releps = 0 +) + +solve_c_parametric( + hypotheses, + test_corr, + alpha, + maxpts = 25000, + abseps = 1e-06, + releps = 0 +) +} +\arguments{ +\item{x}{The root to solve for with \code{stats::uniroot()}.} + +\item{hypotheses}{A numeric vector of hypothesis weights. Must be a vector of +values between 0 & 1 (inclusive). The length should match the length of +\code{p}. The sum of hypothesis weights should not exceed 1.} + +\item{test_corr}{(Optional) A numeric matrix of correlations between test +statistics, which is needed to perform parametric tests using +\code{\link[=adjust_weights_parametric]{adjust_weights_parametric()}}. The number of rows and columns of +this correlation matrix should match the length of \code{p}.} + +\item{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.} + +\item{maxpts}{(Optional) An integer scalar for the maximum number of function +values, which is needed to perform parametric tests using the +\code{mvtnorm::GenzBretz} algorithm. The default is 25000.} + +\item{abseps}{(Optional) A numeric scalar for the absolute error tolerance, +which is needed to perform parametric tests using the \code{mvtnorm::GenzBretz} +algorithm. The default is 1e-6.} + +\item{releps}{(Optional) A numeric scalar for the relative error tolerance +as double, which is needed to perform parametric tests using the +\code{mvtnorm::GenzBretz} algorithm. The default is 0.} + +\item{p}{(Optional) A numeric vector of p-values (unadjusted, raw), whose +values should be between 0 & 1. The length should match the length of +\code{hypotheses}.} + +\item{test_groups}{(Optional) A list of numeric vectors specifying hypotheses +to test together. Grouping is needed to correctly perform Simes and +parametric tests.} +} +\value{ +\itemize{ +\item \code{c_value_function()} returns the difference between +\eqn{\alpha} and the Type I error of the parametric test with the \eqn{c} +value of \code{x}, adjusted for the correlation between test statistics using +parametric tests based on equation (6) of Xi et al. (2017). +\item \code{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). +} +} +\description{ +An intersection hypothesis can be rejected if its p-values are less than or +equal to their adjusted significance levels, which are their adjusted +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 tests: +\itemize{ +\item Parametric tests for \code{\link[=adjust_weights_parametric]{adjust_weights_parametric()}}, +\itemize{ +\item Note that one-sided tests are required for parametric tests. +} +} +} +\references{ +Xi, D., Glimm, E., Maurer, W., and Bretz, F. (2017). A unified framework +for weighted parametric multiple test procedures. +\emph{Biometrical Journal}, 59(5), 918-931. +} +\seealso{ +\code{\link[=adjust_weights_parametric]{adjust_weights_parametric()}} for adjusted hypothesis weights using +parametric tests. +} +\keyword{internal} diff --git a/man/edge_pairs.Rd b/man/edge_pairs.Rd index ca4afd3..6965451 100644 --- a/man/edge_pairs.Rd +++ b/man/edge_pairs.Rd @@ -17,7 +17,4 @@ if no such pairs are found. For an initial graph, find pairs of hypotheses that are connected in both directions. This is used to plot graphs using \code{\link[=plot.initial_graph]{plot.initial_graph()}}. } -\examples{ -graphicalMCP:::edge_pairs(bonferroni_holm(hypotheses = rep(1 / 3, 3))) -} \keyword{internal} diff --git a/man/graph_calculate_power.Rd b/man/graph_calculate_power.Rd index cf8a577..a5f75e6 100644 --- a/man/graph_calculate_power.Rd +++ b/man/graph_calculate_power.Rd @@ -156,11 +156,12 @@ success_fns <- list( ) set.seed(1234) # Bonferroni tests +# Reduce the number of simulations to save time for package compilation power_output <- graph_calculate_power( g, alpha, sim_corr = corr, - sim_n = 1e5, + sim_n = 1e2, power_marginal = marginal_power, sim_success = success_fns ) @@ -174,7 +175,7 @@ graph_calculate_power( test_groups = list(1:2, 3:4), test_types = c("parametric", "simes"), test_corr = list(corr1, NA), - sim_n = 1e3, + sim_n = 1e2, sim_success = list( function(.) .[1] || .[2], function(.) .[1] && .[2] diff --git a/man/graph_create.Rd b/man/graph_create.Rd index 529eae7..2bd9983 100644 --- a/man/graph_create.Rd +++ b/man/graph_create.Rd @@ -101,7 +101,7 @@ transitions <- rbind( H3 = c(0, 1, 0, 0), H4 = c(1, 0, 0, 0) ) -# g <- graph_create(hypotheses, transitions) +\dontrun{g <- graph_create(hypotheses, transitions)} # When names are not specified, hypotheses are numbered sequentially as # H1, H2, ... diff --git a/man/graphicalMCP-package.Rd b/man/graphicalMCP-package.Rd index 0661c9f..07f8c69 100644 --- a/man/graphicalMCP-package.Rd +++ b/man/graphicalMCP-package.Rd @@ -8,7 +8,7 @@ \description{ \if{html}{\figure{logo.png}{options: style='float: right' alt='logo' width='120'}} -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. +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. } \seealso{ Useful links: diff --git a/man/plot.initial_graph.Rd b/man/plot.initial_graph.Rd index fc32c7c..6849762 100644 --- a/man/plot.initial_graph.Rd +++ b/man/plot.initial_graph.Rd @@ -61,7 +61,8 @@ Defaults to all 0, since igraph plots tend to have large margins. It is passed directly to \code{\link[graphics:par]{graphics::par()}} (\code{mar}).} } \value{ -NULL, after plotting the initial graph. +An object x of class \code{initial_graph}, after plotting the initial +graph. } \description{ The plot of an \code{initial_graph} translates the \code{hypotheses} into vertices and diff --git a/man/plot.updated_graph.Rd b/man/plot.updated_graph.Rd index 6ee24c8..1d4b8c5 100644 --- a/man/plot.updated_graph.Rd +++ b/man/plot.updated_graph.Rd @@ -43,7 +43,8 @@ passed directly to \code{\link[graphics:par]{graphics::par()}} (\code{mar}).} }} } \value{ -NULL, after plotting the updated graph. +An object x of class \code{updated_graph}, after plotting the updated +graph. } \description{ Plotting an updated graph is a \emph{very} light wrapper around diff --git a/man/print.graph_report.Rd b/man/print.graph_report.Rd index 948b3b0..6aa3ffb 100644 --- a/man/print.graph_report.Rd +++ b/man/print.graph_report.Rd @@ -19,6 +19,10 @@ to to display.} \item{rows}{An integer scalar indicating how many rows of detailed test results to print.} } +\value{ +An object x of class \code{graph_report}, after printing the report of +conducting a graphical multiple comparison procedure. +} \description{ A printed \code{graph_report} displays the initial graph, p-values and significance levels, rejection decisions, and optional detailed test results. diff --git a/man/print.initial_graph.Rd b/man/print.initial_graph.Rd index 65283b3..a641efb 100644 --- a/man/print.initial_graph.Rd +++ b/man/print.initial_graph.Rd @@ -17,7 +17,8 @@ to to display.} \item{indent}{An integer scalar indicating how many spaces to indent results.} } \value{ -NULL, after printing the initial graph. +An object x of class \code{initial_graph}, after printing the initial +graph. } \description{ A printed \code{initial_graph} displays a header stating "Initial graph", diff --git a/man/print.power_report.Rd b/man/print.power_report.Rd index 544973a..168762f 100644 --- a/man/print.power_report.Rd +++ b/man/print.power_report.Rd @@ -19,6 +19,11 @@ to to display.} \item{rows}{An integer scalar indicating how many rows of detailed test results to print.} } +\value{ +An object x of the class \code{power_report}, after printing the report of +conducting power simulations based on a graphical multiple comparison +procedure. +} \description{ A printed \code{power_report} displays the initial graph, testing and simulation options, power outputs, and optional detailed simulations and test results. diff --git a/man/print.updated_graph.Rd b/man/print.updated_graph.Rd index b36f1ae..f171aac 100644 --- a/man/print.updated_graph.Rd +++ b/man/print.updated_graph.Rd @@ -17,7 +17,8 @@ to to display.} \item{indent}{An integer scalar indicating how many spaces to indent results.} } \value{ -NULL, after printing the updated graph. +An object x of the class \code{updated_graph}, after printing the updated +graph. } \description{ A printed \code{updated_graph} displays the initial graph, the (final) updated