Nothing
Common multiple comparison procedures illustrated using graphicalMCP"
In graphicalMCP: Graphical Multiple Comparison Procedures
knitr::opts_chunk$set(
fig.align = "center",
collapse = TRUE,
comment = "#>",
fig.dim = c(3, 3)
)
library(graphicalMCP)
Introduction
In confirmatory clinical trials, regulatory guidelines mandate the strong control of the family-wise error rate at a prespecified level $\alpha$. Many multiple comparison procedures (MCPs) have been proposed for this purpose. The graphical approaches are a general framework that include many common MCPs as special cases. In this vignette, we illustrate how to use graphicalMCP to perform some common MCPs.
To test $m$ hypotheses using a graphical MCP, each hypothesis $H_i$ receives a weight $0\leq w_i\leq 1$ (called hypothesis weight), where $\sum_{i=1}^{m}w_i\leq 1$. From $H_i$ to $H_j$, there could be a directed and weighted edge $0\leq g_{ij}\leq 1$, which means that when $H_i$ is rejected, its hypothesis weight will be propagated (or transferred) to $H_j$ and $g_{ij}$ determines how much of the propagation. We also require $\sum_{j=1}^{m}g_{ij}\leq 1$ and $g_{ii}=0$.
Bonferroni-based procedures
Bonferroni test
A Bonferroni test splits $\alpha$ equally among hypotheses by testing every hypothesis at a significance level of $\alpha$ divided by the number of hypotheses. Thus it rejects a hypothesis if its p-value is less than or equal to its significance level. There is no propagation between any pair of hypothesis.
set.seed(1234)
alpha <- 0.025
m <- 3
bonferroni_graph <- bonferroni(m)
# transitions <- matrix(0, m, m)
# bonferroni_graph <- graph_create(rep(1 / m, m), transitions)
plot(
bonferroni_graph,
layout = igraph::layout_in_circle(
as_igraph(bonferroni_graph),
order = c(2, 1, 3)
),
vertex.size = 70
)
p_values <- runif(m, 0, alpha)
test_results <-
graph_test_shortcut(
bonferroni_graph,
p = p_values,
alpha = alpha
)
test_results$outputs$rejected
Weighted Bonferroni test
A weighted Bonferroni test splits $\alpha$ among hypotheses by testing every hypothesis at a significance level of $w_i\alpha$. Thus it rejects a hypothesis if its p-value is less than or equal to its significance level. When $w_i=w_j$ for all $i,j$, this means an equal split and the test is the Bonferroni test. There is no propagation between any pair of hypothesis.
set.seed(1234)
alpha <- 0.025
hypotheses <- c(0.5, 0.3, 0.2)
weighted_bonferroni_graph <- bonferroni_weighted(hypotheses)
# m <- length(hypotheses)
# transitions <- matrix(0, m, m)
# weighted_bonferroni_graph <- graph_create(hypotheses, transitions)
plot(
weighted_bonferroni_graph,
layout = igraph::layout_in_circle(
as_igraph(bonferroni_graph),
order = c(2, 1, 3)
),
vertex.size = 70
)
p_values <- runif(m, 0, alpha)
test_results <-
graph_test_shortcut(
weighted_bonferroni_graph,
p = p_values,
alpha = alpha
)
test_results$outputs$rejected
Holm Procedure
Holm (or Bonferroni-Holm) procedures improve over Bonferroni tests by allowing propagation [@holm-1979-simple]. In other words, transition weights between hypotheses may not be zero. So it is uniformly more powerful than Bonferroni tests.
set.seed(1234)
alpha <- 0.025
m <- 3
holm_graph <- bonferroni_holm(m)
# transitions <- matrix(1 / (m - 1), m, m)
# diag(transitions) <- 0
# holm_graph <- graph_create(rep(1 / m, m), transitions)
plot(holm_graph,
layout = igraph::layout_in_circle(
as_igraph(holm_graph),
order = c(2, 1, 3)
),
vertex.size = 70
)
p_values <- runif(m, 0, alpha)
test_results <- graph_test_shortcut(holm_graph, p = p_values, alpha = alpha)
test_results$outputs$rejected
Weighted Holm Procedure
Weighted Holm (or weighted Bonferroni-Holm) procedures improve over weighted Bonferroni tests by allowing propagation [@holm-1979-simple]. In other words, transition weights between hypotheses may not be zero. So it is uniformly more powerful than weighted Bonferroni tests.
set.seed(1234)
alpha <- 0.025
hypotheses <- c(0.5, 0.3, 0.2)
weighted_holm_graph <- bonferroni_holm_weighted(hypotheses)
# m <- length(hypotheses)
# transitions <- matrix(1 / (m - 1), m, m)
# diag(transitions) <- 0
# weighted_holm_graph <- graph_create(hypotheses, transitions)
plot(weighted_holm_graph,
layout = igraph::layout_in_circle(
as_igraph(holm_graph),
order = c(2, 1, 3)
),
vertex.size = 70
)
p_values <- runif(m, 0, alpha)
test_results <- graph_test_shortcut(weighted_holm_graph, p = p_values, alpha = alpha)
test_results$outputs$rejected
Fixed sequence procedure
Fixed sequence (or hierarchical) procedures pre-specify an order of testing [@maurer-1995-multiple,@westfall-2001-optimally]. For example, the procedure will test $H_1$ first. If it is rejected, it will test $H_2$; otherwise the testing stops. If $H_2$ is rejected, it will test $H_3$; otherwise the testing stops. For each hypothesis, it will be tested at the full $\alpha$ level, when it can be tested.
set.seed(1234)
alpha <- 0.025
m <- 3
fixed_sequence_graph <- fixed_sequence(m)
# transitions <- rbind(
# c(0, 1, 0),
# c(0, 0, 1),
# c(0, 0, 0)
# )
# fixed_sequence_graph <- graph_create(c(1, 0, 0), transitions)
plot(fixed_sequence_graph, nrow = 1, asp = 0.5, vertex.size = 50)
p_values <- runif(m, 0, alpha)
test_results <-
graph_test_shortcut(
fixed_sequence_graph,
p = p_values,
alpha = alpha
)
test_results$outputs$rejected
Fallback procedure
Fallback procedures have one-way propagation (like fixed sequence procedures) but allow hypotheses to be tested at different significance levels [@wiens-2003-fixed].
set.seed(1234)
alpha <- 0.025
hypotheses <- c(0.5, 0.3, 0.2)
m <- length(hypotheses)
fallback_graph <- fallback(hypotheses)
# transitions <- rbind(
# c(0, 1, 0),
# c(0, 0, 1),
# c(0, 0, 0)
# )
# fallback_graph <- graph_create(hypotheses, transitions)
plot(fallback_graph, nrow = 1, asp = 0.5, vertex.size = 50)
p_values <- runif(m, 0, alpha)
test_results <-
graph_test_shortcut(
fallback_graph,
p = p_values,
alpha = alpha
)
test_results$outputs$rejected
Further they can be improved to allow propagation from later hypotheses to earlier hypotheses, because it is possible that a later hypothesis is rejected before an earlier hypothesis can be rejected. There are two versions of improvement: fallback_improved_1 due to @wiens-2005-fallback and fallback_improved_2 due to @bretz-2009-graphical respectively.
set.seed(1234)
alpha <- 0.025
hypotheses <- c(0.5, 0.3, 0.2)
m <- length(hypotheses)
fallback_improved_1_graph <- fallback_improved_1(hypotheses)
# transitions <- rbind(
# c(0, 1, 0),
# c(0, 0, 1),
# c(hypotheses[seq_len(m - 1)] / sum(hypotheses[seq_len(m - 1)]), 0)
# )
# fallback_improved_1_graph <- graph_create(hypotheses, transitions)
plot_layout <- rbind(
c(0, 0.8),
c(0.8, 0.8),
c(0.4, 0)
)
plot(
fallback_improved_1_graph,
layout = plot_layout,
asp = 0.5,
edge_curves = c(pairs = 0.5),
vertex.size = 70
)
epsilon <- 0.0001
fallback_improved_2_graph <- fallback_improved_2(hypotheses, epsilon)
# transitions <- rbind(
# c(0, 1, 0),
# c(1 - epsilon, 0, epsilon),
# c(1, 0, 0)
# )
# fallback_improved_2_graph <- graph_create(hypotheses, transitions)
plot_layout <- rbind(
c(0, 0.8),
c(0.8, 0.8),
c(0.4, 0)
)
plot(
fallback_improved_2_graph,
layout = plot_layout,
eps = 0.0001,
asp = 0.5,
edge_curves = c(pairs = 0.5),
vertex.size = 70
)
p_values <- runif(m, 0, alpha)
test_results <-
graph_test_shortcut(
fallback_improved_2_graph,
p = p_values,
alpha = alpha
)
test_results$outputs$rejected
Serial gatekeeping procedure
Serial gatekeeping procedures involve ordered multiple families of hypotheses, where all hypotheses of a family of hypotheses must be rejected before proceeding in the test sequence. The example below considers a primary family consisting of two hypotheses $H_1$ and $H_2$ and a secondary family consisting of a single hypothesis $H_3$. In the primary family, the Holm procedure is applied. If both $H_1$ and $H_2$ are rejected, $H_3$ can be tested at level $\alpha$; otherwise $H_3$ cannot be rejected. To allow the conditional propagation to $H_3$, an $\varepsilon$ edge is used from $H_2$ to $H_3$. It has a very small transition weight so that $H_2$ propagates most of its hypothesis weight to $H_1$ (if not already rejected) and retains a small (non-zero) weight for $H_3$ so that if $H_1$ has been rejected, all hypothesis weight of $H_2$ will be propagated to $H_3$. Here $\varepsilon$ is assigned to be 0.0001 and in practice, the value could be adjusted but it should be much smaller than the smallest p-value observed.
set.seed(1234)
alpha <- 0.025
m <- 3
epsilon <- 0.0001
transitions <- rbind(
c(0, 1, 0),
c(1 - epsilon, 0, epsilon),
c(0, 0, 0)
)
serial_gatekeeping_graph <- graph_create(c(0.5, 0.5, 0), transitions)
plot_layout <- rbind(
c(0, 0.8),
c(0.8, 0.8),
c(0.4, 0)
)
plot(
serial_gatekeeping_graph,
layout = plot_layout,
eps = 0.0001,
asp = 0.5,
edge_curves = c(pairs = 0.5),
vertex.size = 70
)
p_values <- runif(m, 0, alpha)
test_results <-
graph_test_shortcut(
serial_gatekeeping_graph,
p = p_values,
alpha = alpha
)
test_results$outputs$rejected
Parallel gatekeeping procedure
Parallel gatekeeping procedures also involve multiple ordered families of hypotheses, where any null hypotheses of a family of hypotheses must be rejected before proceeding in the test sequence [@dmitrienko-2003-gatekeeping]. The example below considers a primary family consisting of two hypotheses $H_1$ and $H_2$ and a secondary family consisting of two hypotheses $H_3$ and $H_4$. In the primary family, the Bonferroni test is applied. If any of $H_1$ and $H_2$ is rejected, $H_3$ and $H_4$ can be tested at level $\alpha/2$ using the Holm procedure; if both $H_1$ and $H_2$ are rejected, $H_3$ and $H_4$ can be tested at level $\alpha$ using the Holm procedure; otherwise $H_3$ and $H_4$ cannot be
rejected.
set.seed(1234)
alpha <- 0.025
m <- 4
transitions <- rbind(
c(0, 0, 0.5, 0.5),
c(0, 0, 0.5, 0.5),
c(0, 0, 0, 1),
c(0, 0, 1, 0)
)
parallel_gatekeeping_graph <- graph_create(c(0.5, 0.5, 0, 0), transitions)
plot(parallel_gatekeeping_graph, vertex.size = 70)
p_values <- runif(m, 0, alpha)
test_results <-
graph_test_shortcut(
parallel_gatekeeping_graph,
p = p_values,
alpha = alpha
)
test_results$outputs$rejected
The above parallel gatekeeping procedure can be improved by adding $\varepsilon$ edges from secondary hypotheses to primary hypotheses, because it is possible that both secondary hypotheses are rejected but there is still a remaining primary hypothesis
not rejected [@bretz-2009-graphical].
set.seed(1234)
alpha <- 0.025
m <- 4
epsilon <- 0.0001
transitions <- rbind(
c(0, 0, 0.5, 0.5),
c(0, 0, 0.5, 0.5),
c(epsilon, 0, 0, 1 - epsilon),
c(0, epsilon, 1 - epsilon, 0)
)
parallel_gatekeeping_improved_graph <-
graph_create(c(0.5, 0.5, 0, 0), transitions)
plot(parallel_gatekeeping_improved_graph, eps = 0.0001, vertex.size = 70)
p_values <- runif(m, 0, alpha)
test_results <-
graph_test_shortcut(
parallel_gatekeeping_improved_graph,
p = p_values,
alpha = alpha
)
test_results$outputs$rejected
Successive procedure
Successive procedures incorporate successive relationships between hypotheses. For example, the secondary hypothesis is not tested until the primary hypothesis has been rejected. This is similar to using the fixed sequence procedure as a component of a graph. The example below considers two primary hypotheses $H_1$ and $H_2$ and two secondary hypotheses $H_3$ and $H_4$. Primary hypotheses $H_1$ and $H_2$ receive the equal hypothesis weight of 0.5; secondary hypotheses $H_3$ and $H_4$ receive the hypothesis weight of 0. A secondary hypothesis $H_3 (H_4)$ can be tested only if the corresponding primary hypothesis $H_1 (H_2)$ has been rejected. This represents the successive relationships between $H_1$ and $H_3$, and $H_2$ and $H_4$, respectively [@maurer-2011-multiple]. If both $H_1$ and $H_3$ are rejected, their hypothesis weights are propagated to $H_2$ and $H_4$, and vice versa.
set.seed(1234)
alpha <- 0.025
m <- 4
simple_successive_graph <- simple_successive_1()
# transitions <- rbind(
# c(0, 0, 1, 0),
# c(0, 0, 0, 1),
# c(0, 1, 0, 0),
# c(1, 0, 0, 0)
# )
# simple_successive_graph <- graph_create(c(0.5, 0.5, 0, 0), transitions)
plot(simple_successive_graph, layout = "grid", nrow = 2, vertex.size = 70)
p_values <- runif(m, 0, alpha)
test_results <-
graph_test_shortcut(
simple_successive_graph,
p = p_values,
alpha = alpha
)
test_results$outputs$rejected
The above graph could be generalized to allow propagation between primary hypotheses [@maurer-2011-multiple]. A general successive graph is illustrate below with a variable to determine the propagation between $H_1$ and $H_2$.
set.seed(1234)
alpha <- 0.025
m <- 4
successive_var <- simple_successive_var <- function(gamma) {
graph_create(
c(0.5, 0.5, 0, 0),
rbind(
c(0, gamma, 1 - gamma, 0),
c(gamma, 0, 0, 1 - gamma),
c(0, 1, 0, 0),
c(1, 0, 0, 0)
)
)
}
successive_var_graph <- successive_var(0.5)
plot(successive_var_graph, layout = "grid", nrow = 2, vertex.size = 70)
p_values <- runif(m, 0, alpha)
test_results <-
graph_test_shortcut(
successive_var_graph,
p = p_values,
alpha = alpha
)
test_results$outputs$rejected
Hochberg-based procedures
Hochberg procedure
Hochberg procedure [@hochberg-1988-sharper] is a closed test procedure which uses Hochberg tests for every intersection hypothesis. According to @xi-2019-symmetric, the graph for Hochberg procedures is the same as the graph for Holm procedures. Thus to perform Hochberg procedure, we just need to specify test_type to be hochberg.
set.seed(1234)
alpha <- 0.025
m <- 3
hochberg_graph <- hochberg(m)
plot(
hochberg_graph,
layout = igraph::layout_in_circle(
as_igraph(hochberg_graph),
order = c(2, 1, 3)
),
vertex.size = 70
)
p_values <- runif(m, 0, alpha)
test_results <- graph_test_closure(
hochberg_graph,
p = p_values,
alpha = alpha,
test_types = "hochberg"
)
test_results$outputs$rejected
Simes-based procedures
Hommel procedure
Hommel procedure [@hommel-1988-stagewise] is a closed test procedure which uses Simes tests for every intersection hypothesis. According to @xi-2019-symmetric, the graph for Hommel procedures is the same as the graph for Holm procedures. Thus to perform Hommel procedure, we just need to specify test_type to be simes.
set.seed(1234)
alpha <- 0.025
m <- 3
hommel_graph <- hommel(m)
plot(
hommel_graph,
layout = igraph::layout_in_circle(
as_igraph(hommel_graph),
order = c(2, 1, 3)
),
vertex.size = 70
)
p_values <- runif(m, 0, alpha)
test_results <- graph_test_closure(
hommel_graph,
p = p_values,
alpha = alpha,
test_types = "simes"
)
test_results$outputs$rejected
Parametric procedures
Šidák test
The Šidák test is similar to the equally weighted Bonferroni test but it assumes test statistics are independent of each other [@vsidak-1967-rectangular]. Thus it can be performed as a parametric procedure with the identity correlation matrix. Its graph is the same as the graph for the equally weighted Bonferroni test.
set.seed(1234)
alpha <- 0.025
m <- 3
sidak_graph <- sidak(m)
plot(
sidak_graph,
layout = igraph::layout_in_circle(
as_igraph(bonferroni_graph),
order = c(2, 1, 3)
),
vertex.size = 70
)
p_values <- runif(m, 0, alpha)
corr <- diag(m)
test_results <- graph_test_closure(
sidak_graph,
p = p_values, alpha = alpha,
test_types = "parametric",
test_corr = list(corr)
)
test_results$outputs$rejected
Dunnett test
Single-step Dunnett tests are an improvement from Bonferroni tests by incorporating the correlation structure between test statistics [@dunnett-1955-multiple]. Thus their graphs are the same as Bonferroni tests. Assume an equi-correlated case, where the correlation between any pair of test statistics is the same, e.g., 0.5. Then we can perform the single-step Dunnett test by specifying test_type to be parametric and providing the correlation matrix.
set.seed(1234)
alpha <- 0.025
m <- 3
dunnett_graph <- dunnett_single_step(m)
plot(
dunnett_graph,
layout = igraph::layout_in_circle(
as_igraph(dunnett_graph),
order = c(2, 1, 3)
),
vertex.size = 70
)
p_values <- runif(m, 0, alpha)
corr <- matrix(0.5, m, m)
diag(corr) <- 1
test_results <- graph_test_closure(
dunnett_graph,
p = p_values, alpha = alpha,
test_types = "parametric",
test_corr = list(corr)
)
test_results$outputs$rejected
Weighted Dunnett test
Weighted single-step Dunnett tests are an improvement from weighted Bonferroni tests by incorporating the correlation structure between test statistics [@xi-2017-unified]. Thus their graphs are the same as weighted Bonferroni tests. Assume an equi-correlated case, where the correlation between any pair of test statistics is the same, e.g., 0.5. Then we can perform the weighted single-step Dunnett test by specifying test_type to be parametric and providing the correlation matrix.
set.seed(1234)
alpha <- 0.025
hypotheses <- c(0.5, 0.3, 0.2)
dunnett_graph <- dunnett_single_step_weighted(hypotheses)
plot(
dunnett_graph,
layout = igraph::layout_in_circle(
as_igraph(dunnett_graph),
order = c(2, 1, 3)
),
vertex.size = 70
)
p_values <- runif(m, 0, alpha)
corr <- matrix(0.5, m, m)
diag(corr) <- 1
test_results <- graph_test_closure(
dunnett_graph,
p = p_values, alpha = alpha,
test_types = "parametric",
test_corr = list(corr)
)
test_results$outputs$rejected
Dunnett procedure
Dunnett procedures are a closed test procedure and an improvement from Holm procedures by incorporating the correlation structure between test statistics [@xi-2017-unified]. Thus their graphs are the same as Holm procedures. Assume an equi-correlated case, where the correlation between any pair of test statistics is the same, e.g., 0.5. Then we can perform the step-down Dunnett procedure by specifying test_type to be parametric and providing the correlation matrix.
set.seed(1234)
alpha <- 0.025
hypotheses <- c(0.5, 0.3, 0.2)
dunnett_graph <- dunnett_closure_weighted(hypotheses)
plot(
dunnett_graph,
layout = igraph::layout_in_circle(
as_igraph(dunnett_graph),
order = c(2, 1, 3)
),
vertex.size = 70
)
p_values <- runif(m, 0, alpha)
corr <- matrix(0.5, m, m)
diag(corr) <- 1
test_results <- graph_test_closure(
dunnett_graph,
p = p_values, alpha = alpha,
test_types = "parametric",
test_corr = list(corr)
)
test_results$outputs$rejected
Reference
Try the graphicalMCP package in your browser
Any scripts or data that you put into this service are public.
graphicalMCP documentation built on June 8, 2025, 11:19 a.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
knitr::opts_chunk$set( fig.align = "center", collapse = TRUE, comment = "#>", fig.dim = c(3, 3) )
library(graphicalMCP)
Introduction
In confirmatory clinical trials, regulatory guidelines mandate the strong control of the family-wise error rate at a prespecified level $\alpha$. Many multiple comparison procedures (MCPs) have been proposed for this purpose. The graphical approaches are a general framework that include many common MCPs as special cases. In this vignette, we illustrate how to use graphicalMCP to perform some common MCPs.
To test $m$ hypotheses using a graphical MCP, each hypothesis $H_i$ receives a weight $0\leq w_i\leq 1$ (called hypothesis weight), where $\sum_{i=1}^{m}w_i\leq 1$. From $H_i$ to $H_j$, there could be a directed and weighted edge $0\leq g_{ij}\leq 1$, which means that when $H_i$ is rejected, its hypothesis weight will be propagated (or transferred) to $H_j$ and $g_{ij}$ determines how much of the propagation. We also require $\sum_{j=1}^{m}g_{ij}\leq 1$ and $g_{ii}=0$.
Bonferroni-based procedures
Bonferroni test
A Bonferroni test splits $\alpha$ equally among hypotheses by testing every hypothesis at a significance level of $\alpha$ divided by the number of hypotheses. Thus it rejects a hypothesis if its p-value is less than or equal to its significance level. There is no propagation between any pair of hypothesis.
set.seed(1234) alpha <- 0.025 m <- 3 bonferroni_graph <- bonferroni(m) # transitions <- matrix(0, m, m) # bonferroni_graph <- graph_create(rep(1 / m, m), transitions) plot( bonferroni_graph, layout = igraph::layout_in_circle( as_igraph(bonferroni_graph), order = c(2, 1, 3) ), vertex.size = 70 ) p_values <- runif(m, 0, alpha) test_results <- graph_test_shortcut( bonferroni_graph, p = p_values, alpha = alpha ) test_results$outputs$rejected
Weighted Bonferroni test
A weighted Bonferroni test splits $\alpha$ among hypotheses by testing every hypothesis at a significance level of $w_i\alpha$. Thus it rejects a hypothesis if its p-value is less than or equal to its significance level. When $w_i=w_j$ for all $i,j$, this means an equal split and the test is the Bonferroni test. There is no propagation between any pair of hypothesis.
set.seed(1234) alpha <- 0.025 hypotheses <- c(0.5, 0.3, 0.2) weighted_bonferroni_graph <- bonferroni_weighted(hypotheses) # m <- length(hypotheses) # transitions <- matrix(0, m, m) # weighted_bonferroni_graph <- graph_create(hypotheses, transitions) plot( weighted_bonferroni_graph, layout = igraph::layout_in_circle( as_igraph(bonferroni_graph), order = c(2, 1, 3) ), vertex.size = 70 ) p_values <- runif(m, 0, alpha) test_results <- graph_test_shortcut( weighted_bonferroni_graph, p = p_values, alpha = alpha ) test_results$outputs$rejected
Holm Procedure
Holm (or Bonferroni-Holm) procedures improve over Bonferroni tests by allowing propagation [@holm-1979-simple]. In other words, transition weights between hypotheses may not be zero. So it is uniformly more powerful than Bonferroni tests.
set.seed(1234) alpha <- 0.025 m <- 3 holm_graph <- bonferroni_holm(m) # transitions <- matrix(1 / (m - 1), m, m) # diag(transitions) <- 0 # holm_graph <- graph_create(rep(1 / m, m), transitions) plot(holm_graph, layout = igraph::layout_in_circle( as_igraph(holm_graph), order = c(2, 1, 3) ), vertex.size = 70 ) p_values <- runif(m, 0, alpha) test_results <- graph_test_shortcut(holm_graph, p = p_values, alpha = alpha) test_results$outputs$rejected
Weighted Holm Procedure
Weighted Holm (or weighted Bonferroni-Holm) procedures improve over weighted Bonferroni tests by allowing propagation [@holm-1979-simple]. In other words, transition weights between hypotheses may not be zero. So it is uniformly more powerful than weighted Bonferroni tests.
set.seed(1234) alpha <- 0.025 hypotheses <- c(0.5, 0.3, 0.2) weighted_holm_graph <- bonferroni_holm_weighted(hypotheses) # m <- length(hypotheses) # transitions <- matrix(1 / (m - 1), m, m) # diag(transitions) <- 0 # weighted_holm_graph <- graph_create(hypotheses, transitions) plot(weighted_holm_graph, layout = igraph::layout_in_circle( as_igraph(holm_graph), order = c(2, 1, 3) ), vertex.size = 70 ) p_values <- runif(m, 0, alpha) test_results <- graph_test_shortcut(weighted_holm_graph, p = p_values, alpha = alpha) test_results$outputs$rejected
Fixed sequence procedure
Fixed sequence (or hierarchical) procedures pre-specify an order of testing [@maurer-1995-multiple,@westfall-2001-optimally]. For example, the procedure will test $H_1$ first. If it is rejected, it will test $H_2$; otherwise the testing stops. If $H_2$ is rejected, it will test $H_3$; otherwise the testing stops. For each hypothesis, it will be tested at the full $\alpha$ level, when it can be tested.
set.seed(1234) alpha <- 0.025 m <- 3 fixed_sequence_graph <- fixed_sequence(m) # transitions <- rbind( # c(0, 1, 0), # c(0, 0, 1), # c(0, 0, 0) # ) # fixed_sequence_graph <- graph_create(c(1, 0, 0), transitions) plot(fixed_sequence_graph, nrow = 1, asp = 0.5, vertex.size = 50) p_values <- runif(m, 0, alpha) test_results <- graph_test_shortcut( fixed_sequence_graph, p = p_values, alpha = alpha ) test_results$outputs$rejected
Fallback procedure
Fallback procedures have one-way propagation (like fixed sequence procedures) but allow hypotheses to be tested at different significance levels [@wiens-2003-fixed].
set.seed(1234) alpha <- 0.025 hypotheses <- c(0.5, 0.3, 0.2) m <- length(hypotheses) fallback_graph <- fallback(hypotheses) # transitions <- rbind( # c(0, 1, 0), # c(0, 0, 1), # c(0, 0, 0) # ) # fallback_graph <- graph_create(hypotheses, transitions) plot(fallback_graph, nrow = 1, asp = 0.5, vertex.size = 50) p_values <- runif(m, 0, alpha) test_results <- graph_test_shortcut( fallback_graph, p = p_values, alpha = alpha ) test_results$outputs$rejected
Further they can be improved to allow propagation from later hypotheses to earlier hypotheses, because it is possible that a later hypothesis is rejected before an earlier hypothesis can be rejected. There are two versions of improvement: fallback_improved_1 due to @wiens-2005-fallback and fallback_improved_2 due to @bretz-2009-graphical respectively.
set.seed(1234) alpha <- 0.025 hypotheses <- c(0.5, 0.3, 0.2) m <- length(hypotheses) fallback_improved_1_graph <- fallback_improved_1(hypotheses) # transitions <- rbind( # c(0, 1, 0), # c(0, 0, 1), # c(hypotheses[seq_len(m - 1)] / sum(hypotheses[seq_len(m - 1)]), 0) # ) # fallback_improved_1_graph <- graph_create(hypotheses, transitions) plot_layout <- rbind( c(0, 0.8), c(0.8, 0.8), c(0.4, 0) ) plot( fallback_improved_1_graph, layout = plot_layout, asp = 0.5, edge_curves = c(pairs = 0.5), vertex.size = 70 ) epsilon <- 0.0001 fallback_improved_2_graph <- fallback_improved_2(hypotheses, epsilon) # transitions <- rbind( # c(0, 1, 0), # c(1 - epsilon, 0, epsilon), # c(1, 0, 0) # ) # fallback_improved_2_graph <- graph_create(hypotheses, transitions) plot_layout <- rbind( c(0, 0.8), c(0.8, 0.8), c(0.4, 0) ) plot( fallback_improved_2_graph, layout = plot_layout, eps = 0.0001, asp = 0.5, edge_curves = c(pairs = 0.5), vertex.size = 70 ) p_values <- runif(m, 0, alpha) test_results <- graph_test_shortcut( fallback_improved_2_graph, p = p_values, alpha = alpha ) test_results$outputs$rejected
Serial gatekeeping procedure
Serial gatekeeping procedures involve ordered multiple families of hypotheses, where all hypotheses of a family of hypotheses must be rejected before proceeding in the test sequence. The example below considers a primary family consisting of two hypotheses $H_1$ and $H_2$ and a secondary family consisting of a single hypothesis $H_3$. In the primary family, the Holm procedure is applied. If both $H_1$ and $H_2$ are rejected, $H_3$ can be tested at level $\alpha$; otherwise $H_3$ cannot be rejected. To allow the conditional propagation to $H_3$, an $\varepsilon$ edge is used from $H_2$ to $H_3$. It has a very small transition weight so that $H_2$ propagates most of its hypothesis weight to $H_1$ (if not already rejected) and retains a small (non-zero) weight for $H_3$ so that if $H_1$ has been rejected, all hypothesis weight of $H_2$ will be propagated to $H_3$. Here $\varepsilon$ is assigned to be 0.0001 and in practice, the value could be adjusted but it should be much smaller than the smallest p-value observed.
set.seed(1234) alpha <- 0.025 m <- 3 epsilon <- 0.0001 transitions <- rbind( c(0, 1, 0), c(1 - epsilon, 0, epsilon), c(0, 0, 0) ) serial_gatekeeping_graph <- graph_create(c(0.5, 0.5, 0), transitions) plot_layout <- rbind( c(0, 0.8), c(0.8, 0.8), c(0.4, 0) ) plot( serial_gatekeeping_graph, layout = plot_layout, eps = 0.0001, asp = 0.5, edge_curves = c(pairs = 0.5), vertex.size = 70 ) p_values <- runif(m, 0, alpha) test_results <- graph_test_shortcut( serial_gatekeeping_graph, p = p_values, alpha = alpha ) test_results$outputs$rejected
Parallel gatekeeping procedure
Parallel gatekeeping procedures also involve multiple ordered families of hypotheses, where any null hypotheses of a family of hypotheses must be rejected before proceeding in the test sequence [@dmitrienko-2003-gatekeeping]. The example below considers a primary family consisting of two hypotheses $H_1$ and $H_2$ and a secondary family consisting of two hypotheses $H_3$ and $H_4$. In the primary family, the Bonferroni test is applied. If any of $H_1$ and $H_2$ is rejected, $H_3$ and $H_4$ can be tested at level $\alpha/2$ using the Holm procedure; if both $H_1$ and $H_2$ are rejected, $H_3$ and $H_4$ can be tested at level $\alpha$ using the Holm procedure; otherwise $H_3$ and $H_4$ cannot be rejected.
set.seed(1234) alpha <- 0.025 m <- 4 transitions <- rbind( c(0, 0, 0.5, 0.5), c(0, 0, 0.5, 0.5), c(0, 0, 0, 1), c(0, 0, 1, 0) ) parallel_gatekeeping_graph <- graph_create(c(0.5, 0.5, 0, 0), transitions) plot(parallel_gatekeeping_graph, vertex.size = 70) p_values <- runif(m, 0, alpha) test_results <- graph_test_shortcut( parallel_gatekeeping_graph, p = p_values, alpha = alpha ) test_results$outputs$rejected
The above parallel gatekeeping procedure can be improved by adding $\varepsilon$ edges from secondary hypotheses to primary hypotheses, because it is possible that both secondary hypotheses are rejected but there is still a remaining primary hypothesis not rejected [@bretz-2009-graphical].
set.seed(1234) alpha <- 0.025 m <- 4 epsilon <- 0.0001 transitions <- rbind( c(0, 0, 0.5, 0.5), c(0, 0, 0.5, 0.5), c(epsilon, 0, 0, 1 - epsilon), c(0, epsilon, 1 - epsilon, 0) ) parallel_gatekeeping_improved_graph <- graph_create(c(0.5, 0.5, 0, 0), transitions) plot(parallel_gatekeeping_improved_graph, eps = 0.0001, vertex.size = 70) p_values <- runif(m, 0, alpha) test_results <- graph_test_shortcut( parallel_gatekeeping_improved_graph, p = p_values, alpha = alpha ) test_results$outputs$rejected
Successive procedure
Successive procedures incorporate successive relationships between hypotheses. For example, the secondary hypothesis is not tested until the primary hypothesis has been rejected. This is similar to using the fixed sequence procedure as a component of a graph. The example below considers two primary hypotheses $H_1$ and $H_2$ and two secondary hypotheses $H_3$ and $H_4$. Primary hypotheses $H_1$ and $H_2$ receive the equal hypothesis weight of 0.5; secondary hypotheses $H_3$ and $H_4$ receive the hypothesis weight of 0. A secondary hypothesis $H_3 (H_4)$ can be tested only if the corresponding primary hypothesis $H_1 (H_2)$ has been rejected. This represents the successive relationships between $H_1$ and $H_3$, and $H_2$ and $H_4$, respectively [@maurer-2011-multiple]. If both $H_1$ and $H_3$ are rejected, their hypothesis weights are propagated to $H_2$ and $H_4$, and vice versa.
set.seed(1234) alpha <- 0.025 m <- 4 simple_successive_graph <- simple_successive_1() # transitions <- rbind( # c(0, 0, 1, 0), # c(0, 0, 0, 1), # c(0, 1, 0, 0), # c(1, 0, 0, 0) # ) # simple_successive_graph <- graph_create(c(0.5, 0.5, 0, 0), transitions) plot(simple_successive_graph, layout = "grid", nrow = 2, vertex.size = 70) p_values <- runif(m, 0, alpha) test_results <- graph_test_shortcut( simple_successive_graph, p = p_values, alpha = alpha ) test_results$outputs$rejected
The above graph could be generalized to allow propagation between primary hypotheses [@maurer-2011-multiple]. A general successive graph is illustrate below with a variable to determine the propagation between $H_1$ and $H_2$.
set.seed(1234) alpha <- 0.025 m <- 4 successive_var <- simple_successive_var <- function(gamma) { graph_create( c(0.5, 0.5, 0, 0), rbind( c(0, gamma, 1 - gamma, 0), c(gamma, 0, 0, 1 - gamma), c(0, 1, 0, 0), c(1, 0, 0, 0) ) ) } successive_var_graph <- successive_var(0.5) plot(successive_var_graph, layout = "grid", nrow = 2, vertex.size = 70) p_values <- runif(m, 0, alpha) test_results <- graph_test_shortcut( successive_var_graph, p = p_values, alpha = alpha ) test_results$outputs$rejected
Hochberg-based procedures
Hochberg procedure
Hochberg procedure [@hochberg-1988-sharper] is a closed test procedure which uses Hochberg tests for every intersection hypothesis. According to @xi-2019-symmetric, the graph for Hochberg procedures is the same as the graph for Holm procedures. Thus to perform Hochberg procedure, we just need to specify test_type to be hochberg.
set.seed(1234) alpha <- 0.025 m <- 3 hochberg_graph <- hochberg(m) plot( hochberg_graph, layout = igraph::layout_in_circle( as_igraph(hochberg_graph), order = c(2, 1, 3) ), vertex.size = 70 ) p_values <- runif(m, 0, alpha) test_results <- graph_test_closure( hochberg_graph, p = p_values, alpha = alpha, test_types = "hochberg" ) test_results$outputs$rejected
Simes-based procedures
Hommel procedure
Hommel procedure [@hommel-1988-stagewise] is a closed test procedure which uses Simes tests for every intersection hypothesis. According to @xi-2019-symmetric, the graph for Hommel procedures is the same as the graph for Holm procedures. Thus to perform Hommel procedure, we just need to specify test_type to be simes.
set.seed(1234) alpha <- 0.025 m <- 3 hommel_graph <- hommel(m) plot( hommel_graph, layout = igraph::layout_in_circle( as_igraph(hommel_graph), order = c(2, 1, 3) ), vertex.size = 70 ) p_values <- runif(m, 0, alpha) test_results <- graph_test_closure( hommel_graph, p = p_values, alpha = alpha, test_types = "simes" ) test_results$outputs$rejected
Parametric procedures
Šidák test
The Šidák test is similar to the equally weighted Bonferroni test but it assumes test statistics are independent of each other [@vsidak-1967-rectangular]. Thus it can be performed as a parametric procedure with the identity correlation matrix. Its graph is the same as the graph for the equally weighted Bonferroni test.
set.seed(1234) alpha <- 0.025 m <- 3 sidak_graph <- sidak(m) plot( sidak_graph, layout = igraph::layout_in_circle( as_igraph(bonferroni_graph), order = c(2, 1, 3) ), vertex.size = 70 ) p_values <- runif(m, 0, alpha) corr <- diag(m) test_results <- graph_test_closure( sidak_graph, p = p_values, alpha = alpha, test_types = "parametric", test_corr = list(corr) ) test_results$outputs$rejected
Dunnett test
Single-step Dunnett tests are an improvement from Bonferroni tests by incorporating the correlation structure between test statistics [@dunnett-1955-multiple]. Thus their graphs are the same as Bonferroni tests. Assume an equi-correlated case, where the correlation between any pair of test statistics is the same, e.g., 0.5. Then we can perform the single-step Dunnett test by specifying test_type to be parametric and providing the correlation matrix.
set.seed(1234) alpha <- 0.025 m <- 3 dunnett_graph <- dunnett_single_step(m) plot( dunnett_graph, layout = igraph::layout_in_circle( as_igraph(dunnett_graph), order = c(2, 1, 3) ), vertex.size = 70 ) p_values <- runif(m, 0, alpha) corr <- matrix(0.5, m, m) diag(corr) <- 1 test_results <- graph_test_closure( dunnett_graph, p = p_values, alpha = alpha, test_types = "parametric", test_corr = list(corr) ) test_results$outputs$rejected
Weighted Dunnett test
Weighted single-step Dunnett tests are an improvement from weighted Bonferroni tests by incorporating the correlation structure between test statistics [@xi-2017-unified]. Thus their graphs are the same as weighted Bonferroni tests. Assume an equi-correlated case, where the correlation between any pair of test statistics is the same, e.g., 0.5. Then we can perform the weighted single-step Dunnett test by specifying test_type to be parametric and providing the correlation matrix.
set.seed(1234) alpha <- 0.025 hypotheses <- c(0.5, 0.3, 0.2) dunnett_graph <- dunnett_single_step_weighted(hypotheses) plot( dunnett_graph, layout = igraph::layout_in_circle( as_igraph(dunnett_graph), order = c(2, 1, 3) ), vertex.size = 70 ) p_values <- runif(m, 0, alpha) corr <- matrix(0.5, m, m) diag(corr) <- 1 test_results <- graph_test_closure( dunnett_graph, p = p_values, alpha = alpha, test_types = "parametric", test_corr = list(corr) ) test_results$outputs$rejected
Dunnett procedure
Dunnett procedures are a closed test procedure and an improvement from Holm procedures by incorporating the correlation structure between test statistics [@xi-2017-unified]. Thus their graphs are the same as Holm procedures. Assume an equi-correlated case, where the correlation between any pair of test statistics is the same, e.g., 0.5. Then we can perform the step-down Dunnett procedure by specifying test_type to be parametric and providing the correlation matrix.
set.seed(1234) alpha <- 0.025 hypotheses <- c(0.5, 0.3, 0.2) dunnett_graph <- dunnett_closure_weighted(hypotheses) plot( dunnett_graph, layout = igraph::layout_in_circle( as_igraph(dunnett_graph), order = c(2, 1, 3) ), vertex.size = 70 ) p_values <- runif(m, 0, alpha) corr <- matrix(0.5, m, m) diag(corr) <- 1 test_results <- graph_test_closure( dunnett_graph, p = p_values, alpha = alpha, test_types = "parametric", test_corr = list(corr) ) test_results$outputs$rejected
Reference
Try the graphicalMCP package in your browser
Any scripts or data that you put into this service are public.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.