diff --git a/DESCRIPTION b/DESCRIPTION index 272fa68..273bd4e 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,7 +1,7 @@ Type: Package Package: graphicalMCP Title: Graphical Multiple Comparison Procedures -Version: 0.1.0 +Version: 0.1.1 Authors@R: c( person("Dong", "Xi", , "dong.xi1@gilead.com", role = c("aut", "cre")), person("Ethan", "Brockmann", , "ethan.brockmann@atorusresearch.com", role = "aut"), diff --git a/NEWS.md b/NEWS.md index 03773f7..0e4a3e3 100644 --- a/NEWS.md +++ b/NEWS.md @@ -1,3 +1,7 @@ # graphicalMCP 0.1.0 * First release + +# graphicalMCP 0.1.1 + +* De-duplicate test values column names (#75) diff --git a/R/graph_rejection_orderings.R b/R/graph_rejection_orderings.R index bf74299..9df4195 100644 --- a/R/graph_rejection_orderings.R +++ b/R/graph_rejection_orderings.R @@ -44,7 +44,7 @@ graph_rejection_orderings <- function(shortcut_test_result) { rev(expand.grid(rep(list(rejected), length(rejected)))), 1, function(row) { - if (length(unique(row)) == length(row)){ + if (length(unique(row)) == length(row)) { structure(row, names = hyp_names[row]) } else { NULL @@ -61,7 +61,6 @@ graph_rejection_orderings <- function(shortcut_test_result) { intermediate_graph <- graph for (hyp_num in hyp_ordering) { - if (p[[hyp_num]] <= intermediate_graph$hypotheses[[hyp_num]] * alpha) { intermediate_graph <- graph_update(intermediate_graph, hyp_num)$updated_graph diff --git a/R/graph_test_closure.R b/R/graph_test_closure.R index 42d2432..1ef5f6e 100644 --- a/R/graph_test_closure.R +++ b/R/graph_test_closure.R @@ -314,7 +314,8 @@ graph_test_closure <- function(graph, # "c" value is only used in parametric testing, so there's no need to # include this column when there are no parametric groups if (!any(test_types == "parametric")) { - df_test_values[c("c_value", "*")] <- NULL + df_test_values[c("c_value", "*_a")] <- NULL + names(df_test_values)[names(df_test_values) == "*_b"] <- "*" } } diff --git a/R/graph_test_shortcut.R b/R/graph_test_shortcut.R index e6e4e88..f3517ee 100644 --- a/R/graph_test_shortcut.R +++ b/R/graph_test_shortcut.R @@ -146,7 +146,7 @@ graph_test_shortcut <- function(graph, # algorithm") names(test_values_step)[[1]] <- "Step" test_values_step$Step <- step_num - test_values_step[c("Test", "c_value", "*")] <- NULL + test_values_step[c("Test", "c_value", "*_a")] <- NULL df_test_values <- rbind(df_test_values, test_values_step) diff --git a/R/plot.initial_graph.R b/R/plot.initial_graph.R index 8e21d34..f879364 100644 --- a/R/plot.initial_graph.R +++ b/R/plot.initial_graph.R @@ -64,12 +64,12 @@ #' #' epsilon <- 1e-5 #' transitions <- rbind( -#' c(0, 0.5, 0.25, 0, 0.25, 0), -#' c(0.5, 0, 0, 0.25, 0, 0.25), -#' c(0, 0, 0, 0, 1, 0), -#' c(epsilon, 0, 0, 0, 0, 1 - epsilon), -#' c(0, epsilon, 1 - epsilon, 0, 0, 0), -#' c(0, 0, 0, 1, 0, 0) +#' c(0, 0.5, 0.25, 0, 0.25, 0), +#' c(0.5, 0, 0, 0.25, 0, 0.25), +#' c(0, 0, 0, 0, 1, 0), +#' c(epsilon, 0, 0, 0, 0, 1 - epsilon), +#' c(0, epsilon, 1 - epsilon, 0, 0, 0), +#' c(0, 0, 0, 1, 0, 0) #' ) #' #' g <- graph_create(hypotheses, transitions) @@ -77,10 +77,10 @@ #' plot_layout <- rbind( #' c(.15, .5), #' c(.65, .5), -#' c( 0, 0), -#' c( .5, 0), -#' c( .3, 0), -#' c( .8, 0) +#' c(0, 0), +#' c(.5, 0), +#' c(.3, 0), +#' c(.8, 0) #' ) #' #' plot(g, layout = plot_layout, eps = epsilon, edge_curves = c(pairs = .5)) diff --git a/R/print.graph_report.R b/R/print.graph_report.R index 88b4711..bb00aeb 100644 --- a/R/print.graph_report.R +++ b/R/print.graph_report.R @@ -199,6 +199,8 @@ print.graph_report <- function(x, ..., precision = 4, indent = 2, rows = 10) { format(as.numeric(num_col), digits = precision) } ) + + names(crit_res)[grep("\\*_", names(crit_res))] <- "*" if (any(x$inputs$test_types == "parametric")) { crit_res$c_value <- ifelse( trimws(crit_res$c_value) == "NA", diff --git a/R/test_values.R b/R/test_values.R index 45fa91a..a0cb1a0 100644 --- a/R/test_values.R +++ b/R/test_values.R @@ -37,9 +37,9 @@ test_values_bonferroni <- function(p, hypotheses, alpha, intersection = NA) { p = p, "<=" = "<=", c_value = "", - "*" = "", + "*_a" = "", Weight = hypotheses, - "*" = "*", + "*_b" = "*", Alpha = alpha, Inequality_holds = ifelse( p == 0 & hypotheses == 0, @@ -69,9 +69,9 @@ test_values_parametric <- function(p, p = p, "<=" = "<=", c_value = c_value, - "*" = "*", + "*_a" = "*", Weight = hypotheses, - "*" = "*", + "*_b" = "*", Alpha = alpha, Inequality_holds = ifelse( p == 0 & hypotheses == 0, @@ -103,9 +103,9 @@ test_values_simes <- function(p, hypotheses, alpha, intersection = NA) { p = p, "<=" = "<=", c_value = "", - "*" = "", + "*_a" = "", Weight = w_sum, - "*" = "*", + "*_b" = "*", Alpha = alpha, Inequality_holds = ifelse( p == 0 & w_sum == 0, diff --git a/man/graph_calculate_power.Rd b/man/graph_calculate_power.Rd index 0a1cf4e..ce6ec2d 100644 --- a/man/graph_calculate_power.Rd +++ b/man/graph_calculate_power.Rd @@ -57,7 +57,7 @@ A list with three elements \item power - A list of measures of how often hypotheses are rejected \itemize{ \item power_local - Rejection proportion for each hypothesis individually -\item power_expected - Average number of hypotheses rejected in a single +\item rejection_expected - Average number of hypotheses rejected in a single simulation \item power_at_least_1 - Proportion of simulations which reject any hypothesis diff --git a/man/graphicalMCP-package.Rd b/man/graphicalMCP-package.Rd index 4b8117c..0661c9f 100644 --- a/man/graphicalMCP-package.Rd +++ b/man/graphicalMCP-package.Rd @@ -4,11 +4,11 @@ \name{graphicalMCP-package} \alias{graphicalMCP} \alias{graphicalMCP-package} -\title{graphicalMCP: Graphical Approach for Multiple Comparison Procedures} +\title{graphicalMCP: Graphical Multiple Comparison Procedures} \description{ \if{html}{\figure{logo.png}{options: style='float: right' alt='logo' width='120'}} -A multiple comparison procedure (or multiple test procedure) is a a statistical analysis method for determining efficacy of multiple drugs, or multiple doses of the same drug, in a single clinical trial. (Bretz et al., 2011) laid out a graph-based approach to these multiple comparison procedures, in which the weights of the vertices and edges of the graph are determined independently of the particular statistical test used to assess results. This is a low-dependency implementation of the methods described in that and subsequent papers. +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. } \seealso{ Useful links: diff --git a/tests/testthat/test-graph_test_closure.R b/tests/testthat/test-graph_test_closure.R index 863e2b9..e1b45f9 100644 --- a/tests/testthat/test-graph_test_closure.R +++ b/tests/testthat/test-graph_test_closure.R @@ -308,7 +308,6 @@ test_that("compare adjusted p-values to lrstat - Bonferroni & Simes", { ignore_attr = TRUE ) } - } }) diff --git a/vignettes/closed-testing.Rmd b/vignettes/closed-testing.Rmd index 68f662c..0476ed0 100644 --- a/vignettes/closed-testing.Rmd +++ b/vignettes/closed-testing.Rmd @@ -80,12 +80,12 @@ Given a set of p-values for $H_1, \ldots, H_6$, the graphical multiple compariso ```{r graph-test-bonferroni} p_values <- c(0.015, 0.013, 0.01, 0.007, 0.1, 0.0124) test_results <- graph_test_closure( - g, - p = p_values, - alpha = 0.025, - verbose = TRUE, + g, + p = p_values, + alpha = 0.025, + verbose = TRUE, test_values = TRUE - ) +) test_results$outputs$adjusted_p @@ -203,7 +203,7 @@ sample_size <- rep(200, 3) unpooled_variance <- prop[-1] * (1 - prop[-1]) / sample_size[-1] + - prop[1] * (1 - prop[1]) / sample_size[1] + prop[1] * (1 - prop[1]) / sample_size[1] noncentrality_parameter_primary <- -(prop[-1] - prop[1]) / sqrt(unpooled_variance) @@ -230,7 +230,7 @@ variance <- sd[-1]^2 / sample_size[-1] + sd[1]^2 / sample_size[1] noncentrality_parameter_se1 <- (mean_change_se1[-1] - mean_change_se1[1]) / - sqrt(variance) + sqrt(variance) marginal_power_se1 <- pnorm( qnorm(alpha, lower.tail = FALSE), @@ -246,7 +246,7 @@ mean_change_se2 <- c(6, 8, 9) noncentrality_parameter_se2 <- (mean_change_se2[-1] - mean_change_se2[1]) / - sqrt(variance) + sqrt(variance) marginal_power_se2 <- pnorm( qnorm(alpha, lower.tail = FALSE), @@ -309,20 +309,20 @@ These are the default outputs from the `graph_calculate_power` function. In addi success_fns <- list( # Probability to reject H1 H1 = function(x) x[1], - + # Expected number of rejections `Expected no. of rejections` = function(x) x[1] + x[2] + x[3] + x[4] + x[5] + x[6], # Probability to reject at least one hypothesis `AtLeast1` = function(x) x[1] | x[2] | x[3] | x[4] | x[5] | x[6], - + # Probability to reject all hypotheses `All` = function(x) x[1] & x[2] & x[3] & x[4] & x[5] & x[6], - + # Probability to reject both H1 and H2 `H1andH2` = function(x) x[1] & x[2], - + # Probability to reject all hypotheses for the low dose or the high dose `(H1andH3andH5)or(H2andH4andH6)` <- function(x) (x[1] & x[3] & x[5]) | (x[2] & x[4] & x[6]) diff --git a/vignettes/comparisons.Rmd b/vignettes/comparisons.Rmd index c9f4264..86490eb 100644 --- a/vignettes/comparisons.Rmd +++ b/vignettes/comparisons.Rmd @@ -60,11 +60,11 @@ for (i in 1:1000) { graph <- random_graph(5) graphicalmcp_weights <- graphicalMCP::graph_generate_weights(graph) dimnames(graphicalmcp_weights) <- list(NULL, NULL) - gmcp_weights <- + gmcp_weights <- gMCP::generateWeights(graph$transitions, graph$hypotheses) gmcp_weights <- gmcp_weights[nrow(gmcp_weights):1, ] # Reorder rows identical <- c( - identical, + identical, all.equal(gmcp_weights, graphicalmcp_weights, tolerance = 1e-7) ) } @@ -83,12 +83,12 @@ identical <- NULL for (i in 1:10000) { graph <- random_graph(5) p <- runif(5, 0, alpha) - graphicalmcp_test_shortcut <- + graphicalmcp_test_shortcut <- graph_test_shortcut(graph, p, alpha = alpha)$outputs$adjusted_p - gmcp_test_shortcut <- + gmcp_test_shortcut <- gMCP(as_graphMCP(graph), p, alpha = alpha)@adjPValues identical <- c( - identical, + identical, all.equal(graphicalmcp_test_shortcut, gmcp_test_shortcut, tolerance = 1e-7) ) } @@ -110,9 +110,9 @@ for (i in 1:1000) { graph <- random_graph(5) marginal_power <- runif(5, 0.5, 0.9) noncentrality_parameter <- - qnorm(1 - 0.025, lower.tail = TRUE) - - qnorm(1 - marginal_power, lower.tail = TRUE) - + qnorm(1 - 0.025, lower.tail = TRUE) - + qnorm(1 - marginal_power, lower.tail = TRUE) + set.seed(1234 + i - 1) graphicalmcp_power <- rbind( graphicalmcp_power, @@ -124,7 +124,7 @@ for (i in 1:1000) { sim_n = 2^17 )$power$power_local ) - + set.seed(1234 + i - 1) gmcp_power <- rbind( gmcp_power, @@ -143,7 +143,7 @@ diff <- data.frame( rbind(graphicalmcp_power, gmcp_power), procedure = rep(c("graphicalMCP", "gMCP"), each = nrow(graphicalmcp_power)) ) - + write.csv( diff, here::here("vignettes/cache/comparisons_power_shortcut.csv"), @@ -156,10 +156,10 @@ gmcp_power <- subset(diff, procedure == "gMCP") round( max( abs( - graphicalmcp_power[, - ncol(graphicalmcp_power)] - - gmcp_power[, - ncol(gmcp_power)] - ) - ), + graphicalmcp_power[, -ncol(graphicalmcp_power)] - + gmcp_power[, -ncol(gmcp_power)] + ) + ), 4 ) # Maximum difference in local power among 1000 cases ``` @@ -191,8 +191,8 @@ test_corr <- rbind( for (i in 1:10000) { p <- runif(4, 0, alpha) graphicalmcp_test_parametric <- graph_test_closure( - graph, - p, + graph, + p, alpha = alpha, test_groups = list(1:2, 3:4), test_types = c("parametric", "bonferroni"), @@ -205,7 +205,7 @@ for (i in 1:10000) { correlation = test_corr )@adjPValues identical <- c( - identical, + identical, all.equal(graphicalmcp_test_parametric, gmcp_test_parametric, tolerance = 1e-7) ) } @@ -241,7 +241,7 @@ for (i in 1:100) { marginal_power <- runif(4, 0.5, 0.9) noncentrality_parameter <- qnorm(1 - 0.025, lower.tail = TRUE) - - qnorm(1 - marginal_power, lower.tail = TRUE) + qnorm(1 - marginal_power, lower.tail = TRUE) set.seed(1234 + i - 1) graphicalmcp_power_parametric <- rbind( @@ -257,7 +257,7 @@ for (i in 1:100) { sim_n = 2^14 )$power$power_local ) - + set.seed(1234 + i - 1) gmcp_power_parametric <- rbind( gmcp_power_parametric, @@ -277,7 +277,7 @@ diff <- data.frame( rbind(graphicalmcp_power_parametric, gmcp_power_parametric), procedure = rep(c("graphicalMCP", "gMCP"), each = nrow(gmcp_power_parametric)) ) - + write.csv( diff, here::here("vignettes/cache/comparisons_power_parametric.csv"), @@ -290,10 +290,10 @@ gmcp_power <- subset(diff, procedure == "gMCP") round( max( abs( - graphicalmcp_power_parametric[, - ncol(graphicalmcp_power_parametric)] - - gmcp_power_parametric[, - ncol(gmcp_power)] - ) - ), 4 + graphicalmcp_power_parametric[, -ncol(graphicalmcp_power_parametric)] - + gmcp_power_parametric[, -ncol(gmcp_power)] + ) + ), 4 ) # Maximum difference in local power among 100 cases ``` @@ -324,8 +324,8 @@ family <- rbind( for (i in 1:10000) { p <- runif(4, 0, alpha) graphicalmcp_test_simes <- graph_test_closure( - graph, - p, + graph, + p, alpha = alpha, test_groups = list(1:2, 3:4), test_types = c("simes", "bonferroni") @@ -336,10 +336,10 @@ for (i in 1:10000) { fwgtmat(graph$hypotheses, graph$transitions), p, family - ) - + ) + identical <- c( - identical, + identical, all.equal(graphicalmcp_test_simes, lrstat_test_simes, tolerance = 1e-7) ) } diff --git a/vignettes/generate-closure.Rmd b/vignettes/generate-closure.Rmd index 90e5e7f..78eb3e9 100644 --- a/vignettes/generate-closure.Rmd +++ b/vignettes/generate-closure.Rmd @@ -218,10 +218,10 @@ gcp_conventional <- function( num_hyps <- length(graph$hypotheses) graphicalMCP:::power_input_val( - graph, - sim_n, - power_marginal, - sim_corr, + graph, + sim_n, + power_marginal, + sim_corr, sim_success ) @@ -361,29 +361,29 @@ benchmark_list <- lapply( # A graph for the Holm procedure transitions <- matrix(1 / (num_hyps - 1), num_hyps, num_hyps) diag(transitions) <- 0 - + graph <- graph_create(rep(1 / num_hyps, num_hyps), transitions) - + corr <- diag(num_hyps) - + benchmark <- mark( `gMCP (C)` = gMCP::calcPower( - graph$hypotheses, - .025, + graph$hypotheses, + .025, graph$transitions, - n.sim = 2^14, + n.sim = 2^14, corr.sim = corr ), - `graphicalMCP conventional (R)` = + `graphicalMCP conventional (R)` = gcp_conventional( - graph, + graph, .025, power_marginal = rep(.9, num_hyps), sim_n = 2^14 ), `graphicalMCP parent-child (R)` = graph_calculate_power( - graph, + graph, .025, power_marginal = rep(.9, num_hyps), sim_n = 2^14 @@ -393,11 +393,11 @@ benchmark_list <- lapply( time_unit = "s", min_iterations = 5 )[, 1:5] - + benchmark <- benchmark[benchmark$median > .00001, ] benchmark$char_expression <- as.character(benchmark$expression) benchmark <- benchmark[, c("char_expression", "median")] - + cbind(data.frame(num_hyps = num_hyps), benchmark) } ) @@ -430,25 +430,25 @@ benchmark_list <- lapply( graph <- fixed_sequence(num_hyps) corr <- diag(num_hyps) - + benchmark <- mark( `gMCP (C)` = gMCP::calcPower( - graph$hypotheses, - .025, + graph$hypotheses, + .025, graph$transitions, - n.sim = 2^14, + n.sim = 2^14, corr.sim = corr ), - `graphicalMCP conventional (R)` = + `graphicalMCP conventional (R)` = gcp_conventional( - graph, + graph, .025, power_marginal = rep(.9, num_hyps), sim_n = 2^14 ), `graphicalMCP parent-child (R)` = graph_calculate_power( - graph, + graph, .025, power_marginal = rep(.9, num_hyps), sim_n = 2^14 @@ -458,11 +458,11 @@ benchmark_list <- lapply( time_unit = "s", min_iterations = 5 )[, 1:5] - + benchmark <- benchmark[benchmark$median > .00001, ] benchmark$char_expression <- as.character(benchmark$expression) benchmark <- benchmark[, c("char_expression", "median")] - + cbind(data.frame(num_hyps = num_hyps), benchmark) } ) @@ -486,7 +486,7 @@ write.csv( ``` ```{r plot-power-benchmarks, fig.dim=c(8, 5), eval=FALSE} -benchmarks_holm <- +benchmarks_holm <- read.csv(here::here("vignettes/cache/power_benchmarks_holm.csv")) benchmarks_fixed_sequence <- read.csv(here::here("vignettes/cache/power_benchmarks_fixed_sequence.csv")) @@ -513,7 +513,7 @@ benchmarks_plot_standard <- geom_line(size = 1) + # scale_y_continuous(labels = scales::label_comma(suffix = "ms")) + scale_color_discrete() + - facet_wrap(~ Procedure, labeller = label_both)+ + facet_wrap(~Procedure, labeller = label_both) + labs( title = "Computing time of power simulations", subtitle = "Median runtime in seconds", diff --git a/vignettes/graphicalMCP.Rmd b/vignettes/graphicalMCP.Rmd index e031160..e81107d 100644 --- a/vignettes/graphicalMCP.Rmd +++ b/vignettes/graphicalMCP.Rmd @@ -72,7 +72,8 @@ Given the set of p-values of all hypotheses, graphical MCPs can be performed usi test_results <- graph_test_shortcut( example_graph, p = c(.01, .005, .03, .01), - alpha = .025) + alpha = .025 +) test_results$outputs$rejected ``` diff --git a/vignettes/shortcut-testing.Rmd b/vignettes/shortcut-testing.Rmd index 7432f6c..8b7cbc6 100644 --- a/vignettes/shortcut-testing.Rmd +++ b/vignettes/shortcut-testing.Rmd @@ -93,7 +93,7 @@ The order of rejections may not be unique and not all orders are valid. For this ```{r graph-test-order} # Obtain all valid orders of rejections -orders <- graph_rejection_orderings(test_results)$valid_orderings +orders <- graph_rejection_orderings(test_results)$valid_orderings orders # Intermediate graphs following the order of H2 and H4 @@ -130,7 +130,7 @@ sample_size <- rep(200, 3) unpooled_variance <- prop[-1] * (1 - prop[-1]) / sample_size[-1] + - prop[1] * (1 - prop[1] ) / sample_size[1] + prop[1] * (1 - prop[1]) / sample_size[1] noncentrality_parameter_primary <- -(prop[-1] - prop[1]) / sqrt(unpooled_variance) @@ -208,19 +208,19 @@ These are the default outputs from the `graph_calculate_power` function. In addi success_fns <- list( # Probability to reject H1 H1 = function(x) x[1], - + # Expected number of rejections `Expected no. of rejections` = function(x) x[1] + x[2] + x[3] + x[4], # Probability to reject at least one hypothesis `AtLeast1` = function(x) x[1] | x[2] | x[3] | x[4], - + # Probability to reject all hypotheses `All` = function(x) x[1] & x[2] & x[3] & x[4], - + # Probability to reject both H1 and H2 `H1andH2` = function(x) x[1] & x[2], - + # Probability to reject both hypotheses for the low dose or the high dose `(H1andH3)or(H2andH4)` = function(x) (x[1] & x[3]) | (x[2] & x[4]) )