Calculate power for a graph-based multiple test procedure using p-values
Source:R/calc_power.R
calc_power_pvals.RdCalculate multiplicity-adjusted (marginal) power values for a graph-based
multiple test procedure applied to a matrix of p-values. Provide both the
hyp_weight vector and transition matrix trans_matrix for the graph you
wish to evaluate.
Usage
calc_power_pvals(
pvals,
hyp_weight,
trans_matrix,
...,
alpha = 0.025,
custom_power = NULL,
sum_to_one_constraint = TRUE,
call = rlang::caller_env()
)Arguments
- pvals
A numeric matrix of p-values. Each row represents a simulated trial, and each column corresponds to a hypothesis.
- hyp_weight
A numeric vector representing the weights of the nodes (hypotheses) in the graph. All elements must be in the interval [0, 1] and the vector must sum to 1.
- trans_matrix
A numeric matrix representing the transition matrix between nodes. All elements must be in the interval [0, 1], with diagonal elements equal to 0, and each row must sum to 1.
- ...
These dots are for future extensions and must be empty.
- alpha
A single numeric value representing the overall one-sided significance level. Default is 0.025.
- custom_power
A list of user-defined power functions. Alternatively, a single function or
multigrain_trial_successobject can be provided. There are two ways to specify this:Anonymous functions can be provided to specify the success criteria. Functions must take one simulation's logical vector of results as an input, and return a scalar. For example, if one is interested in the power to reject hypotheses 1 and 3 one could specify:
f = function(x) {x[1] && x[3]}. If the power of rejecting hypotheses 1 and 2 is also of interest one would use an (optionally named) list:f = list( power1and3 = function(x) {x[1] && x[3]}, power1and2 = function(x) {x[1] && x[2]} ). If the list has no names, the functions will be referenced as"func1","func2", etc. in the output. The user can also provide amultigrain_trial_successobject instead (resulting in a faster calculation).Instead of anonymous functions, one can pass a
multigrain_trial_successobject (or within a list as with the anonymous functions). Seetrial_success()and the examples.
- sum_to_one_constraint
A logical value controlling whether to allow graphs where transition matrix rows are not constrained to sum to one, for example in a fixed sequence. Defaults to
TRUE.- call
The execution environment of a currently running function, e.g.
caller_env(). The function will be mentioned in error messages as the source of the error. See thecallargument ofabort()for more information.
Value
A list containing:
local_power: The local power for each hypothesis: the proportion of simulations in which each hypothesis is rejected.exp_rejections: The expected number of rejections across all simulations.disj_power: Disjunctive power: the probability of rejecting at least one hypothesis.conj_power: Conjunctive power: the probability of rejecting all hypotheses."...": Additional results as specified by user-defined success functions incustom_power.
Details
calc_power_pvals() calculates several power metrics:
local power (probability to reject each individual hypothesis),
disjunctive power (the probability to reject at least one hypothesis),
conjunctive power (the probability to reject all hypotheses), and
the expected number of rejections. Optionally, you can specify user-defined trial success criteria via
custom_power.
Examples
# First we simulate our p-value distribution
power_nominal <- c(0.90, 0.87, 0.73)
corr_matrix <- matrix(
c(
1, 0.2, 0.2,
0.2, 1, 0.2,
0.2, 0.2, 1
),
nrow = 3,
byrow = TRUE
)
pvals <- simulate_pvalues(
power_nominal = power_nominal,
corr_matrix = corr_matrix
)
# Second we construct our graph (in this case a fixed sequence)
hyp_weights <- c(1, 0, 0)
trans_matrix <- matrix(
c(
0, 1, 0,
0, 0, 1,
0, 0, 0
),
nrow = 3,
byrow = TRUE
)
# Third, we construct a list of metrics we wish to evaluate our graph with
trial_success_measure <- trial_success((r1 && r2) || r3)
#> ✔ Trial success function compiled and sourced successfully.
power_metrics <- list(
trial_success = trial_success_measure,
average_power = function(x) {
x[1] + x[2] + x[3]
}
)
# Finally we calculate the power of graph conditional on p-value distribution
result <- calc_power_pvals(
pvals = pvals,
hyp_weight = hyp_weights,
trans_matrix = trans_matrix,
custom_power = power_metrics,
sum_to_one_constraint = FALSE # As third row of graph does not sum to 1
)
result
#> $local_power
#> [1] 0.90018 0.79120 0.60123
#>
#> $exp_rejections
#> [1] 2.29261
#>
#> $disj_power
#> [1] 0.90018
#>
#> $conj_power
#> [1] 0.60123
#>
#> $trial_success
#> [1] 0.7912
#>
#> $average_power
#> [1] 2.29261
#>