Skip to content

Commit

Permalink
Merge branch test structure and graph_update into vignettes-eb
Browse files Browse the repository at this point in the history
This update covers a few structural changes that have been finalized.
The arguments relating to testing strategy now work with names or order,
and `graph_update()` gets a custom deletion order feature.
  • Loading branch information
EeethB committed Sep 6, 2023
2 parents 5315955 + 8bb0d2c commit 153ec21
Show file tree
Hide file tree
Showing 62 changed files with 854 additions and 562 deletions.
2 changes: 1 addition & 1 deletion .github/workflows/pkgdown.yaml
Original file line number Diff line number Diff line change
Expand Up @@ -4,7 +4,7 @@ on:
push:
branches: [vignettes-eb, feature-pkgdown, dev, main, release]
pull_request:
branches: [dev, main, release]
branches: [feature-power, dev, main, release]
release:
types: [published]
workflow_dispatch:
Expand Down
16 changes: 8 additions & 8 deletions R/adjust_p.R
Original file line number Diff line number Diff line change
@@ -1,9 +1,9 @@
#' Calculate adjusted p-values
#'
#' @param p A named numeric vector of p-values
#' @param hypotheses A named numeric vector of hypothesis weights
#' @param corr (Optional) A numeric matrix of correlations between hypotheses'
#' test statistics
#' @param p A numeric vector of p-values
#' @param hypotheses A numeric vector of hypothesis weights
#' @param test_corr (Optional) A numeric matrix of correlations between
#' hypotheses' test statistics
#'
#' @return A single adjusted p-value for the given group
#'
Expand Down Expand Up @@ -45,7 +45,7 @@ adjust_p_bonferroni <- function(p, hypotheses) {
}

#' @rdname adjust_p
adjust_p_parametric <- function(p, hypotheses, corr = NULL) {
adjust_p_parametric <- function(p, hypotheses, test_corr = NULL) {
if (sum(hypotheses) == 0) {
return(Inf)
}
Expand All @@ -60,7 +60,7 @@ adjust_p_parametric <- function(p, hypotheses, corr = NULL) {
1 - mvtnorm::pmvnorm(
lower = -Inf,
upper = z,
corr = corr[w_nonzero, w_nonzero, drop = FALSE],
corr = test_corr[w_nonzero, w_nonzero, drop = FALSE],
algorithm = mvtnorm::GenzBretz(maxpts = 25000, abseps = 1e-6, releps = 0)
)[[1]]
)
Expand Down Expand Up @@ -93,8 +93,8 @@ adjust_p_simes <- function(p, hypotheses) {
# smaller, incorrect adjusted weight (larger, incorrect adjusted p-value).
# The hypothesis that comes second will be correct. [adjust_weights_simes()]
# is only used in power calculations where it should not be possible to have
# identical p-values, since they are sampled randomly (unless `all(corr ==
# 1)`). Furthermore, even when there are incorrect adjusted weights, it
# identical p-values, since they are sampled randomly (unless `all(test_corr
# == 1)`). Furthermore, even when there are incorrect adjusted weights, it
# cannot affect the hypothesis rejections. See Bonferroni function above for
# na.rm reasoning
adjusted_p <- min(
Expand Down
74 changes: 41 additions & 33 deletions R/adjust_weights.R
Original file line number Diff line number Diff line change
Expand Up @@ -2,21 +2,23 @@
#'
#' The weights created by [graph_generate_weights()] work immediately for
#' Bonferroni testing, but parametric and Simes testing require additional
#' calculations. The `adjust_weights_*()` functions apply parametric or
#' Simes weight increases to get updated weights for testing. They also subset
#' the weights columns by the appropriate groups
#' calculations. The `adjust_weights_*()` functions apply parametric or Simes
#' weight increases to get updated weights for testing. They also subset the
#' weights columns by the appropriate groups
#'
#' @param weighting_strategy For parametric, a compact representation of
#' [graph_generate_weights()] output, where missing hypotheses get a missing
#' value for weights, and h-vectors are dropped. For Simes, just the weights
#' from [graph_generate_weights()] output
#' @param corr A numeric matrix of correlations between hypotheses' test
#' @param matrix_weights The second half of columns from
#' [graph_generate_weights()] output, indicating the weights of each
#' intersection
#' @param matrix_intersections The first half of columns from
#' [graph_generate_weights()] output, indicating which hypotheses are
#' contained in each intersection
#' @param test_corr A numeric matrix of correlations between hypotheses' test
#' statistics
#' @param alpha A numeric scalar specifying the global significance level for
#' testing
#' @param p A numeric vector of p-values
#' @param groups A list of numeric vectors specifying hypotheses to test
#' @param test_groups A list of numeric vectors specifying hypotheses to test
#' together
#' @param p A numeric vector of p-values
#' @param hypotheses A numeric vector of hypothesis weights
#' @param x The root to solve for with [stats::uniroot()]
#'
Expand All @@ -42,28 +44,34 @@
#' gw_0 <- gw_large[, 7:12]
#' gw <- ifelse(gw_large[, 1:6], gw_0, NA)
#'
#' graphicalMCP:::adjust_weights_parametric(gw, diag(6), .05, list(1:3))
#' graphicalMCP:::adjust_weights_parametric(
#' gw_0,
#' gw_large[, 1:6],
#' diag(6),
#' .05,
#' list(1:3)
#' )
#'
#' graphicalMCP:::adjust_weights_simes(gw_0, p, list(4:6))
adjust_weights_parametric <- function(weighting_strategy,
corr,
adjust_weights_parametric <- function(matrix_weights,
matrix_intersections,
test_corr,
alpha,
groups) {
matrix_intersections <- !is.na(weighting_strategy)

test_groups) {
c_values <- matrix(
nrow = nrow(weighting_strategy),
ncol = ncol(weighting_strategy),
dimnames = dimnames(weighting_strategy)
nrow = nrow(matrix_weights),
ncol = ncol(matrix_weights),
dimnames = dimnames(matrix_weights)
)

for (group in groups) {
for (row in seq_len(nrow(weighting_strategy))) {
for (group in test_groups) {
for (row in seq_len(nrow(matrix_weights))) {
group_by_intersection <-
group[as.logical(matrix_intersections[row, , drop = TRUE][group])]

group_c_value <- solve_c_parametric(
weighting_strategy[row, group_by_intersection, drop = TRUE],
corr[group_by_intersection, group_by_intersection, drop = FALSE],
matrix_weights[row, group_by_intersection, drop = TRUE],
test_corr[group_by_intersection, group_by_intersection, drop = FALSE],
alpha
)

Expand All @@ -72,22 +80,22 @@ adjust_weights_parametric <- function(weighting_strategy,
}
}

adjusted_weights <- c_values * weighting_strategy
adjusted_weights <- c_values * matrix_weights

adjusted_weights[, unlist(groups), drop = FALSE]
adjusted_weights[, unlist(test_groups), drop = FALSE]
}

#' @rdname adjusted-weights
adjust_weights_simes <- function(weighting_strategy, p, groups) {
adjust_weights_simes <- function(matrix_weights, p, groups) {
ordered_p <- order(p)

weighting_strategy <- weighting_strategy[, ordered_p, drop = FALSE]
matrix_weights <- matrix_weights[, ordered_p, drop = FALSE]

group_adjusted_weights <- vector("list", length(groups))

for (i in seq_along(groups)) {
group_adjusted_weights[[i]] <- matrixStats::rowCumsums(
weighting_strategy[, ordered_p %in% groups[[i]], drop = FALSE],
matrix_weights[, ordered_p %in% groups[[i]], drop = FALSE],
useNames = TRUE
)
}
Expand All @@ -96,7 +104,7 @@ adjust_weights_simes <- function(weighting_strategy, p, groups) {
}

#' @rdname adjusted-weights
c_value_function <- function(x, hypotheses, corr, alpha) {
c_value_function <- function(x, hypotheses, test_corr, alpha) {
hyps_nonzero <- which(hypotheses > 0)
z <- stats::qnorm(x * hypotheses[hyps_nonzero] * alpha, lower.tail = FALSE)
y <- ifelse(
Expand All @@ -105,7 +113,7 @@ c_value_function <- function(x, hypotheses, corr, alpha) {
1 - mvtnorm::pmvnorm(
lower = -Inf,
upper = z,
corr = corr[hyps_nonzero, hyps_nonzero, drop = FALSE],
corr = test_corr[hyps_nonzero, hyps_nonzero, drop = FALSE],
algorithm = mvtnorm::GenzBretz(maxpts = 25000, abseps = 1e-6, releps = 0)
)[[1]]
)
Expand All @@ -114,7 +122,7 @@ c_value_function <- function(x, hypotheses, corr, alpha) {
}

#' @rdname adjusted-weights
solve_c_parametric <- function(hypotheses, corr, alpha) {
solve_c_parametric <- function(hypotheses, test_corr, alpha) {
num_hyps <- seq_along(hypotheses)
c_value <- ifelse(
length(num_hyps) == 1 || sum(hypotheses) == 0,
Expand All @@ -123,12 +131,12 @@ solve_c_parametric <- function(hypotheses, corr, alpha) {
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(corr == 1)
# 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,
corr = corr,
test_corr = test_corr,
alpha = alpha
)$root
)
Expand Down
4 changes: 2 additions & 2 deletions R/example_graphs.R
Original file line number Diff line number Diff line change
Expand Up @@ -2,8 +2,8 @@
#'
#' @param num_hyps Number of vertices in the graph
#' @param hypotheses Hypothesis weights in a fallback procedure
#' @param hyp_names Optional names for the hypotheses (Must have length
#' `num_hyps` or be NULL)
#' @param hyp_names (Optional) A character vector of hypothesis names. If names
#' are not specified, hypotheses will be named sequentially as H1, H2, ...
#'
#' @return An S3 object as returned by [graph_create()]
#'
Expand Down
Loading

0 comments on commit 153ec21

Please sign in to comment.