Update rejection orderings function to be modified graph report

openpharma · Sep 7, 2023 · 0d12012 · 0d12012
1 parent 947a35c
commit 0d12012
Show file tree

Hide file tree

Showing 12 changed files with 91 additions and 31 deletions.
diff --git a/NAMESPACE b/NAMESPACE
@@ -21,6 +21,7 @@ export(fixed_sequence)
 export(graph_calculate_power)
 export(graph_create)
 export(graph_generate_weights)
+export(graph_rejection_orderings)
 export(graph_test_closure)
 export(graph_test_shortcut)
 export(graph_update)

diff --git a/R/graph_rejection_orderings.R b/R/graph_rejection_orderings.R
@@ -1,3 +1,30 @@
+#' Find alternate rejection orderings for shortcut testing
+#'
+#' When using `graph_test_shortcut()`, there may be multiple hypotheses which
+#' can be rejected at each step. The specific hypothesis chosen is decided based
+#' on the minimum adjusted p-value. This function shows every other order that
+#' deletes the same hypotheses, where each hypothesis is still a valid deletion
+#' at the step it is chosen.
+#'
+#' @param shortcut_test_result A `graph_report` object as returned by
+#'   [graph_test_shortcut()]
+#'
+#' @return A modified `graph_report` object containing all valid orderings for
+#'   deleting the significant hypotheses
+#' @export
+#'
+#' @examples
+#' graph <- simple_successive_2()
+#'
+#' short_res <- graph_test_shortcut(graph, c(.018, .01, .03, .004))
+#'
+#' # Reject H1, H2, and H4
+#' short_res$outputs$rejected
+#'
+#' # But these cannot be rejected in any order - For instance, H4 has 0 weight
+#' # in the initial graph and cannot be rejected first
+#'
+#' graph_rejection_orderings(short_res)$valid_orderings
 graph_rejection_orderings <- function(shortcut_test_result) {
   # Extract basic testing values -----------------------------------------------
   graph <- shortcut_test_result$inputs$graph
@@ -7,44 +34,54 @@ graph_rejection_orderings <- function(shortcut_test_result) {
   hyp_names <- names(graph$hypotheses)
 
   # Permute rejected hypotheses ------------------------------------------------
-  rejected_names <- hyp_names[shortcut_test_result$outputs$rejected]
+  rejected <- which(shortcut_test_result$outputs$rejected)
 
   list_possible_orderings <- apply(
-    expand.grid(rep(list(rejected_names), length(rejected_names))),
+    expand.grid(rep(list(rejected), length(rejected))),
     1,
     function(row) if (length(unique(row)) == length(row)) unname(row) else NULL
   )
   list_possible_orderings <- Filter(Negate(is.null), list_possible_orderings)
 
   # Find which permutations are valid rejection orderings ----------------------
   orderings_valid <- vector("logical", length(list_possible_orderings))
-  graph_sequences <- rep(
-    list(c(list(graph), vector("list", length(rejected_names)))),
-    length(list_possible_orderings)
-  )
 
   for (hyp_ordering_num in seq_along(list_possible_orderings)) {
     hyp_ordering <- list_possible_orderings[[hyp_ordering_num]]
     intermediate_graph <- graph
-    graph_index <- 2
 
-    for (hyp_name in hyp_ordering) {
-      hyp_num <- which(hyp_name == names(graph$hypotheses))
+    for (hyp_num in hyp_ordering) {
 
       if (p[[hyp_num]] <= intermediate_graph$hypotheses[[hyp_num]] * alpha) {
         intermediate_graph <-
-          graph_update(intermediate_graph, hyp_name == hyp_names)$updated_graph
-
-        graph_sequences[[hyp_ordering_num]][[graph_index]] <- intermediate_graph
+          graph_update(intermediate_graph, hyp_num)$updated_graph
       } else {
         orderings_valid[[hyp_ordering_num]] <- FALSE
         break
       }
 
       orderings_valid[[hyp_ordering_num]] <- TRUE
-      graph_index <- graph_index + 1
     }
   }
 
-  list(valid_hypothesis_sequences = list_possible_orderings[orderings_valid])
+  list_orderings_code <- lapply(
+    list_possible_orderings[orderings_valid],
+    function(ordering) {
+      paste0("c(", paste(ordering, collapse = ", "), ")")
+    }
+  )
+
+  valid_orderings <- do.call(rbind, list_orderings_code)
+  dimnames(valid_orderings) <- list(
+    seq_len(length(list_orderings_code)),
+    "Valid rejection orderings"
+  )
+
+  structure(
+    c(
+      shortcut_test_result,
+      list(valid_orderings = valid_orderings)
+    ),
+    class = "graph_report"
+  )
 }
diff --git a/R/plot.initial_graph.R b/R/plot.initial_graph.R
@@ -34,7 +34,7 @@
 #' plot(simple_successive_2(), layout = "grid")
 plot.initial_graph <- function(x,
                                ...,
-                               layout = igraph::layout_nicely,
+                               layout = "grid",
                                nrow = NULL,
                                ncol = NULL,
                                edge_curves = NULL,

diff --git a/R/print.graph_report.R b/R/print.graph_report.R
@@ -233,6 +233,17 @@ print.graph_report <- function(x, ..., precision = 4, indent = 2, rows = 10) {
     }
   }
 
+  # Optional alternate orderings
+  if (!is.null(x$valid_orderings)) {
+    section_break("Alternate rejection orderings ($valid_rejection_orderings)")
+
+    apply(
+      x$valid_orderings, 1,
+      function(ordering) cat(pad, ordering, "\n", sep = "")
+    )
+    cat("\n")
+  }
+
   invisible(x)
 }
 

diff --git a/README.Rmd b/README.Rmd
@@ -63,7 +63,7 @@ transitions <- rbind(
 hyp_names <- c("H1", "H2", "H3", "H4")
 example_graph <- graph_create(hypotheses, transitions, hyp_names)
 
-plot(example_graph, layout = "grid", asp = .7)
+plot(example_graph, layout = "grid")
 ```
 
 ### Update graph
@@ -124,11 +124,11 @@ These methods were originally implemented in the [gMCP package](https://github.c
 
 However, because development has ceased on the original package, we hope to re-implement the methods with a more general testing framework; with fewer dependencies, in particular shedding the Java dependency; with simpler, more transparent S3 classes; and with improvements to the accuracy of the parametric and Simes test methods.
 
-A portion of Simes testing is also implemented in the lrstat package (`install.packages("lrstat")`).
+A portion of Simes testing is also implemented in the lrstat package.
 
 ## Citation
 
-```{r citation}
+```{r citation, results="markup"}
 citation("graphicalMCP")
 ```
 
@@ -175,9 +175,9 @@ To that end, there are several entities encountered in the world of graphical MC
 +---------------------------------+--------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------+---------------------------------------+---------------------------------+---------------------------------+
 | **Testing strategy**            | **Test types** and **test groups** combined                                                                                                                                                                        |                                       |                                 |                                 |
 +---------------------------------+--------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------+---------------------------------------+---------------------------------+---------------------------------+
-| **Marginal power**              | The power to reject each null **hypothesis** at full **alpha**. Closely related to the non-centrality parameter, which is the mean of each null **hypothesis** in the underlying multivariate normal distribution: |                                       | `power_marginal`                | Correlation matrix              |
+| **Marginal power**              | The power to reject each null **hypothesis** at full **alpha**. Closely related to the non-centrality parameter, which is the mean of each null **hypothesis** in the underlying multivariate normal distribution: |                                       | `marginal_power`                | Correlation matrix              |
 |                                 |                                                                                                                                                                                                                    |                                       |                                 |                                 |
-|                                 | `ncp = qnorm(1 - alpha) - qnorm(1 - power_marginal)`                                                                                                                                                               |                                       |                                 |                                 |
+|                                 | `ncp = qnorm(1 - alpha) - qnorm(1 - marginal_power)`                                                                                                                                                               |                                       |                                 |                                 |
 +---------------------------------+--------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------+---------------------------------------+---------------------------------+---------------------------------+
 | **Correlation matrix**          | Specification of correlations between **hypotheses**. Together with **marginal power**, this specifies the (known or assumed) underlying multivariate normal distribution of the null **hypotheses**.              |                                       | `corr`, `test_corr`, `sim_corr` | Marginal power                  |
 +---------------------------------+--------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------+---------------------------------------+---------------------------------+---------------------------------+

diff --git a/README.md b/README.md
@@ -66,7 +66,7 @@ transitions <- rbind(
 hyp_names <- c("H1", "H2", "H3", "H4")
 example_graph <- graph_create(hypotheses, transitions, hyp_names)
 
-plot(example_graph, layout = "grid", asp = .7)
+plot(example_graph, layout = "grid")
 ```
 
 <img src="man/figures/README-create-graph-1.png" width="100%" />
@@ -161,13 +161,26 @@ fewer dependencies, in particular shedding the Java dependency; with
 simpler, more transparent S3 classes; and with improvements to the
 accuracy of the parametric and Simes test methods.
 
-A portion of Simes testing is also implemented in the lrstat package
-(`install.packages("lrstat")`).
+A portion of Simes testing is also implemented in the lrstat package.
 
 ## Citation
 
 ``` r
 citation("graphicalMCP")
+#> To cite graphicalMCP in publications use:
+#> 
+#>   Xi, D.; Brockmann, E. (2023). graphicalMCP: Graph-based multiple
+#>   comparison procedures. version 0.1.0. Gilead Sciences, Inc. Foster
+#>   City, California. https://github.com/Gilead-BioStats/graphicalMCP
+#> 
+#>   Frank Bretz, Martin Posch, Ekkehard Glimm, Florian Klinglmueller,
+#>   Willi Maurer, Kornelius Rohmeyer (2011), Graphical approaches for
+#>   multiple comparison procedures using weighted Bonferroni, Simes or
+#>   parametric tests. Biometrical Journal 53 (6), pages 894--913, Wiley.
+#> 
+#> To see these entries in BibTeX format, use 'print(<citation>,
+#> bibtex=TRUE)', 'toBibtex(.)', or set
+#> 'options(citation.bibtex.max=999)'.
 ```
 
 ## Acknowledgments
@@ -332,9 +345,9 @@ combined</td>
 <strong>alpha</strong>. Closely related to the non-centrality parameter,
 which is the mean of each null <strong>hypothesis</strong> in the
 underlying multivariate normal distribution:</p>
-<p><code>ncp = qnorm(1 - alpha) - qnorm(1 - power_marginal)</code></p></td>
+<p><code>ncp = qnorm(1 - alpha) - qnorm(1 - marginal_power)</code></p></td>
 <td></td>
-<td><code>power_marginal</code></td>
+<td><code>marginal_power</code></td>
 <td>Correlation matrix</td>
 </tr>
 <tr class="even">

diff --git a/_pkgdown.yml b/_pkgdown.yml
@@ -18,6 +18,7 @@ reference:
 - contents:
   - graph_update
   - print.updated_graph
+  - plot.updated_graph
 - title: Calculating sub-graph weights
 - desc: Functions for generating the weights of the closure
 - contents:

diff --git a/man/figures/README-create-graph-1.png b/man/figures/README-create-graph-1.png
diff --git a/man/plot.initial_graph.Rd b/man/plot.initial_graph.Rd
diff --git a/man/plot.updated_graph.Rd b/man/plot.updated_graph.Rd
diff --git a/misc/plot_igraph_complex.R b/misc/plot_igraph_complex.R
@@ -21,7 +21,6 @@ plot(
   layout = complex_layout,
   edge_curves = c("pairs" = 1),
   vertex.size = 15,
-  asp = .4,
   edge.label.cex = 1.2,
   edge.arrow.size = 1
   )
@@ -50,7 +49,6 @@ plot(
                   "H3|H7" = .05,
                   "H9|H1" = .05),
   vertex.size = 15,
-  asp = .4,
   edge.label.cex = 1.2,
   edge.arrow.size = 1
 )
diff --git a/vignettes/graph-examples.Rmd b/vignettes/graph-examples.Rmd
@@ -47,7 +47,7 @@ fixed_sequence_graph <- graph_create(
   )
 )
 
-plot(fixed_sequence_graph, layout = "grid", nrow = 1, asp = .1)
+plot(fixed_sequence_graph, layout = "grid", nrow = 1)
 ```
 
 # Simple successive
@@ -134,7 +134,6 @@ plot(
   wd_2005_graph,
   layout = "grid",
   nrow = 1,
-  asp = .1,
   edge_curves = c(pairs = -6, "H3|H1" = -6)
 )
 ```