## ----include = FALSE---------------------------------------------------------- knitr::opts_chunk$set( fig.align = "center", collapse = TRUE, comment = "#>", cache.lazy = FALSE ) ## ----setup, message = FALSE, warning = FALSE---------------------------------- library(graphicalMCP) ## ----bonferroni-results------------------------------------------------------- out <- read.csv(here::here("vignettes/internal-validation_bonferroni.csv")) # Matching power using the adjusted p-value approach all.equal(out$adjusted_p, rep(TRUE, nrow(out))) # Matching power using the adjusted significance level approach all.equal(out$adjusted_significance_level, rep(TRUE, nrow(out))) ## ----hochberg-results--------------------------------------------------------- out <- read.csv(here::here("vignettes/internal-validation_hochberg.csv")) # Matching power using the adjusted p-value approach all.equal(out$adjusted_p, rep(TRUE, nrow(out))) # Matching power using the adjusted significance level approach all.equal(out$adjusted_significance_level, rep(TRUE, nrow(out))) ## ----simes-results------------------------------------------------------------ out <- read.csv(here::here("vignettes/internal-validation_simes.csv")) # Matching power using the adjusted p-value approach all.equal(out$adjusted_p, rep(TRUE, nrow(out))) # Matching power using the adjusted significance level approach all.equal(out$adjusted_significance_level, rep(TRUE, nrow(out))) ## ----parametric-results------------------------------------------------------- out <- read.csv(here::here("vignettes/internal-validation_parametric.csv")) # Matching power using the adjusted p-value approach all.equal(out$adjusted_p, rep(TRUE, nrow(out))) # Matching power using the adjusted significance level approach all.equal(out$adjusted_significance_level, rep(TRUE, nrow(out))) ## ----mixed-results------------------------------------------------------------ out <- read.csv(here::here("vignettes/internal-validation_mixed.csv")) # Matching power using the adjusted p-value approach all.equal(out$adjusted_p, rep(TRUE, nrow(out))) # Matching power using the adjusted significance level approach all.equal(out$adjusted_significance_level, rep(TRUE, nrow(out))) ## ----parametric-mixed-results------------------------------------------------- out <- read.csv(here::here("vignettes/internal-validation_parametric-mixed.csv")) # Matching power using the adjusted p-value approach all.equal(out$adjusted_p, rep(TRUE, nrow(out))) # Matching power using the adjusted significance level approach all.equal(out$adjusted_significance_level, rep(TRUE, nrow(out))) ## ----bonferroni, include = FALSE, eval = FALSE-------------------------------- # n_graph <- 1000 # m <- 6 # alpha <- 0.025 # n_sim <- 1e2 # sim_corr <- matrix(0.5, m, m) # diag(sim_corr) <- 1 # out_adjusted_p <- rep(NA, n_graph) # out_adjusted_significance_level <- rep(NA, n_graph) # for (i in 1:n_graph) { # set.seed(1234 + i - 1) # graph <- random_graph(m) # marginal_power <- runif(m, 0.5, 0.9) # results_calculate_power <- graph_calculate_power( # graph, # alpha = alpha, # power_marginal = marginal_power, # sim_corr = sim_corr, # sim_n = n_sim, # verbose = TRUE # ) # p_sim <- results_calculate_power$details$p_sim # adjusted_p_results <- p_sim # adjusted_significance_level_results <- p_sim # for (j in 1:n_sim) { # temp <- graph_test_shortcut( # graph, # p_sim[j, ], # alpha = alpha, # test_values = TRUE # ) # adjusted_p_results[j, ] <- temp$outputs$rejected # temp_adjusted_significance_level <- temp$test_values$results$Inequality_holds # names(temp_adjusted_significance_level) <- temp$test_values$results$Hypothesis # adjusted_significance_level_results[j, ] <- # temp_adjusted_significance_level[names(temp$outputs$rejected)] # } # results_adjusted_p <- colMeans(adjusted_p_results) # names(results_adjusted_p) <- names(temp$outputs$rejected) # # Matching power using the adjusted p-value approach # out_adjusted_p[i] <- all.equal( # results_calculate_power$power$power_local, # results_adjusted_p, # tolerance = 1 / n_sim # ) # # results_significance_level <- colMeans(adjusted_significance_level_results) # names(results_significance_level) <- names(temp$outputs$rejected) # # Matching power using the adjusted significance level approach # out_adjusted_significance_level[i] <- all.equal( # results_calculate_power$power$power_local, # results_significance_level, # tolerance = 1 / n_sim # ) # } # out <- cbind(out_adjusted_p, out_adjusted_significance_level) # colnames(out) <- c("adjusted_p", "adjusted_significance_level") # write.csv( # out, # here::here("vignettes/internal-validation_bonferroni.csv"), # row.names = FALSE # ) ## ----hochberg, include = FALSE, eval = FALSE---------------------------------- # n_graph <- 1000 # m <- 4 # alpha <- 0.025 # n_sim <- 1e2 # sim_corr <- matrix(0.5, m, m) # diag(sim_corr) <- 1 # out_adjusted_p <- rep(NA, n_graph) # out_adjusted_significance_level <- rep(NA, n_graph) # for (i in 1:n_graph) { # set.seed(1234 + i - 1) # groups <- sample(1:m) # test_groups <- list(groups[1:(m / 2)], groups[(m / 2 + 1):m]) # graph <- random_graph(m) # graph$hypotheses[groups[1:(m / 2)]] <- sum(graph$hypotheses[groups[1:(m / 2)]]) / 3 # graph$hypotheses[groups[(m / 2 + 1):m]] <- sum(graph$hypotheses[groups[(m / 2 + 1):m]]) / 3 # graph$transitions <- matrix(1 / (m - 1), nrow = m, ncol = m) # diag(graph$transitions) <- 0 # marginal_power <- runif(m, 0.5, 0.9) # results_calculate_power <- graph_calculate_power( # graph, # alpha = alpha, # power_marginal = marginal_power, # sim_corr = sim_corr, # sim_n = n_sim, # test_types = c("hochberg", "hochberg"), # test_groups = test_groups, # verbose = TRUE # ) # p_sim <- results_calculate_power$details$p_sim # adjusted_p_results <- p_sim # adjusted_significance_level_results <- p_sim # for (j in 1:n_sim) { # temp <- graph_test_closure( # graph, # p_sim[j, ], # alpha = alpha, # test_types = c("hochberg", "hochberg"), # test_groups = test_groups, # test_values = TRUE # ) # adjusted_p_results[j, ] <- temp$outputs$rejected # intersection <- unique(temp$test_values$results$Intersection) # intersection_rej <- intersection # for (k in 1:length(intersection)) { # intersection_rej[k] <- any(subset( # temp$test_values$results, # Intersection == intersection[k] # )$Inequality_holds) # } # for (k in 1:m) { # id <- substr(intersection, start = k, stop = k) == "1" # adjusted_significance_level_results[j, k] <- suppressWarnings(all(intersection_rej[id])) # } # } # results_adjusted_p <- colMeans(adjusted_p_results) # names(results_adjusted_p) <- names(temp$outputs$rejected) # # Matching power using the adjusted p-value approach # out_adjusted_p[i] <- all.equal( # results_calculate_power$power$power_local, # results_adjusted_p, # tolerance = 1 / n_sim # ) # # results_significance_level <- colMeans(adjusted_significance_level_results) # names(results_significance_level) <- names(temp$outputs$rejected) # # Matching power using the adjusted significance level approach # out_adjusted_significance_level[i] <- all.equal( # results_calculate_power$power$power_local, # results_significance_level, # tolerance = 1 / n_sim # ) # } # out <- cbind(out_adjusted_p, out_adjusted_significance_level) # colnames(out) <- c("adjusted_p", "adjusted_significance_level") # write.csv( # out, # here::here("vignettes/internal-validation_hochberg.csv"), # row.names = FALSE # ) ## ----simes, include = FALSE, eval = FALSE------------------------------------- # n_graph <- 1000 # m <- 4 # alpha <- 0.025 # n_sim <- 1e2 # sim_corr <- matrix(0.5, m, m) # diag(sim_corr) <- 1 # out_adjusted_p <- rep(NA, n_graph) # out_adjusted_significance_level <- rep(NA, n_graph) # for (i in 1:n_graph) { # set.seed(1234 + i - 1) # groups <- sample(1:m) # test_groups <- list(groups[1:(m / 2)], groups[(m / 2 + 1):m]) # graph <- random_graph(m) # marginal_power <- runif(m, 0.5, 0.9) # results_calculate_power <- graph_calculate_power( # graph, # alpha = alpha, # power_marginal = marginal_power, # sim_corr = sim_corr, # sim_n = n_sim, # test_types = c("simes", "simes"), # test_groups = test_groups, # verbose = TRUE # ) # p_sim <- results_calculate_power$details$p_sim # adjusted_p_results <- p_sim # adjusted_significance_level_results <- p_sim # for (j in 1:n_sim) { # temp <- graph_test_closure( # graph, # p_sim[j, ], # alpha = alpha, # test_types = c("simes", "simes"), # test_groups = test_groups, # test_values = TRUE # ) # adjusted_p_results[j, ] <- temp$outputs$rejected # intersection <- unique(temp$test_values$results$Intersection) # intersection_rej <- intersection # for (k in 1:length(intersection)) { # intersection_rej[k] <- any(subset( # temp$test_values$results, # Intersection == intersection[k] # )$Inequality_holds) # } # for (k in 1:m) { # id <- substr(intersection, start = k, stop = k) == "1" # adjusted_significance_level_results[j, k] <- suppressWarnings(all(intersection_rej[id])) # } # } # results_adjusted_p <- colMeans(adjusted_p_results) # names(results_adjusted_p) <- names(temp$outputs$rejected) # # Matching power using the adjusted p-value approach # out_adjusted_p[i] <- all.equal( # results_calculate_power$power$power_local, # results_adjusted_p, # tolerance = 1 / n_sim # ) # # results_significance_level <- colMeans(adjusted_significance_level_results) # names(results_significance_level) <- names(temp$outputs$rejected) # # Matching power using the adjusted significance level approach # out_adjusted_significance_level[i] <- all.equal( # results_calculate_power$power$power_local, # results_significance_level, # tolerance = 1 / n_sim # ) # } # out <- cbind(out_adjusted_p, out_adjusted_significance_level) # colnames(out) <- c("adjusted_p", "adjusted_significance_level") # write.csv( # out, # here::here("vignettes/internal-validation_simes.csv"), # row.names = FALSE # ) ## ----parametric, include = FALSE, eval = FALSE-------------------------------- # n_graph <- 1000 # m <- 4 # alpha <- 0.025 # n_sim <- 1e2 # sim_corr <- matrix(0.5, m, m) # diag(sim_corr) <- 1 # out_adjusted_p <- rep(NA, n_graph) # out_adjusted_significance_level <- rep(NA, n_graph) # for (i in 1:n_graph) { # set.seed(1234 + i - 1) # groups <- sample(1:m) # test_groups <- list(groups[1:(m / 2)], groups[(m / 2 + 1):m]) # test_corr <- list( # sim_corr[groups[1:(m / 2)], groups[1:(m / 2)]], # sim_corr[groups[(m / 2 + 1):m], groups[(m / 2 + 1):m]] # ) # graph <- random_graph(m) # marginal_power <- runif(m, 0.5, 0.9) # results_calculate_power <- graph_calculate_power( # graph, # alpha = alpha, # power_marginal = marginal_power, # sim_corr = sim_corr, # sim_n = n_sim, # test_types = c("parametric", "parametric"), # test_groups = test_groups, # test_corr = test_corr, # verbose = TRUE # ) # p_sim <- results_calculate_power$details$p_sim # adjusted_p_results <- p_sim # adjusted_significance_level_results <- p_sim # for (j in 1:n_sim) { # temp <- graph_test_closure( # graph, # p_sim[j, ], # alpha = alpha, # test_types = c("parametric", "parametric"), # test_groups = test_groups, # test_corr = test_corr, # test_values = TRUE # ) # adjusted_p_results[j, ] <- temp$outputs$rejected # intersection <- unique(temp$test_values$results$Intersection) # intersection_rej <- intersection # for (k in 1:length(intersection)) { # intersection_rej[k] <- any(subset( # temp$test_values$results, # Intersection == intersection[k] # )$Inequality_holds) # } # for (k in 1:m) { # id <- substr(intersection, start = k, stop = k) == "1" # adjusted_significance_level_results[j, k] <- suppressWarnings(all(intersection_rej[id])) # } # } # results_adjusted_p <- colMeans(adjusted_p_results) # names(results_adjusted_p) <- names(temp$outputs$rejected) # # Matching power using the adjusted p-value approach # out_adjusted_p[i] <- all.equal( # results_calculate_power$power$power_local, # results_adjusted_p, # tolerance = 1 / n_sim # ) # # results_significance_level <- colMeans(adjusted_significance_level_results) # names(results_significance_level) <- names(temp$outputs$rejected) # # Matching power using the adjusted significance level approach # out_adjusted_significance_level[i] <- all.equal( # results_calculate_power$power$power_local, # results_significance_level, # tolerance = 1 / n_sim # ) # } # out <- cbind(out_adjusted_p, out_adjusted_significance_level) # colnames(out) <- c("adjusted_p", "adjusted_significance_level") # write.csv( # out, # here::here("vignettes/internal-validation_parametric.csv"), # row.names = FALSE # ) ## ----mixed, include = FALSE, eval = FALSE------------------------------------- # n_graph <- 1000 # m <- 4 # alpha <- 0.025 # n_sim <- 1e2 # sim_corr <- matrix(0.5, m, m) # diag(sim_corr) <- 1 # out_adjusted_p <- rep(NA, n_graph) # out_adjusted_significance_level <- rep(NA, n_graph) # for (i in 1:n_graph) { # set.seed(1234 + i - 1) # groups <- sample(1:m) # test_groups <- list(groups[1:(m / 2)], groups[(m / 2 + 1):m]) # test_corr <- list( # sim_corr[groups[1:(m / 2)], groups[1:(m / 2)]], # sim_corr[groups[(m / 2 + 1):m], groups[(m / 2 + 1):m]] # ) # test_types <- sample(c("bonferroni", "hochberg", "simes"), 2) # graph <- random_graph(m) # graph$hypotheses[groups[1:(m / 2)]] <- sum(graph$hypotheses[groups[1:(m / 2)]]) / 3 # graph$hypotheses[groups[(m / 2 + 1):m]] <- sum(graph$hypotheses[groups[(m / 2 + 1):m]]) / 3 # graph$transitions <- matrix(1 / (m - 1), nrow = m, ncol = m) # diag(graph$transitions) <- 0 # marginal_power <- runif(m, 0.5, 0.9) # results_calculate_power <- graph_calculate_power( # graph, # alpha = alpha, # power_marginal = marginal_power, # sim_corr = sim_corr, # sim_n = n_sim, # test_types = test_types, # test_groups = test_groups, # verbose = TRUE # ) # p_sim <- results_calculate_power$details$p_sim # adjusted_p_results <- p_sim # adjusted_significance_level_results <- p_sim # for (j in 1:n_sim) { # temp <- graph_test_closure( # graph, # p_sim[j, ], # alpha = alpha, # test_types = test_types, # test_groups = test_groups, # test_values = TRUE # ) # adjusted_p_results[j, ] <- temp$outputs$rejected # intersection <- unique(temp$test_values$results$Intersection) # intersection_rej <- intersection # for (k in 1:length(intersection)) { # intersection_rej[k] <- any(subset( # temp$test_values$results, # Intersection == intersection[k] # )$Inequality_holds) # } # for (k in 1:m) { # id <- substr(intersection, start = k, stop = k) == "1" # adjusted_significance_level_results[j, k] <- suppressWarnings(all(intersection_rej[id])) # } # } # results_adjusted_p <- colMeans(adjusted_p_results) # names(results_adjusted_p) <- names(temp$outputs$rejected) # # Matching power using the adjusted p-value approach # out_adjusted_p[i] <- all.equal( # results_calculate_power$power$power_local, # results_adjusted_p, # tolerance = 1 / n_sim # ) # # results_significance_level <- colMeans(adjusted_significance_level_results) # names(results_significance_level) <- names(temp$outputs$rejected) # # Matching power using the adjusted significance level approach # out_adjusted_significance_level[i] <- all.equal( # results_calculate_power$power$power_local, # results_significance_level, # tolerance = 1 / n_sim # ) # } # out <- cbind(out_adjusted_p, out_adjusted_significance_level) # colnames(out) <- c("adjusted_p", "adjusted_significance_level") # write.csv( # out, # here::here("vignettes/internal-validation_mixed.csv"), # row.names = FALSE # ) ## ----parametric-mixed, include = FALSE, eval = FALSE-------------------------- # n_graph <- 1000 # m <- 4 # alpha <- 0.025 # n_sim <- 1e2 # sim_corr <- matrix(0.5, m, m) # diag(sim_corr) <- 1 # out_adjusted_p <- rep(NA, n_graph) # out_adjusted_significance_level <- rep(NA, n_graph) # for (i in 1:n_graph) { # set.seed(1234 + i - 1) # groups <- sample(1:m) # test_groups <- list(groups[1:(m / 2)], groups[(m / 2 + 1):m]) # test_corr <- list(sim_corr[groups[1:(m / 2)], groups[1:(m / 2)]], NA) # test_types <- c( # "parametric", # sample(c("bonferroni", "hochberg", "simes"), 1) # ) # graph <- random_graph(m) # graph$hypotheses[groups[1:(m / 2)]] <- sum(graph$hypotheses[groups[1:(m / 2)]]) / 3 # graph$hypotheses[groups[(m / 2 + 1):m]] <- sum(graph$hypotheses[groups[(m / 2 + 1):m]]) / 3 # graph$transitions <- matrix(1 / (m - 1), nrow = m, ncol = m) # diag(graph$transitions) <- 0 # marginal_power <- runif(m, 0.5, 0.9) # results_calculate_power <- graph_calculate_power( # graph, # alpha = alpha, # power_marginal = marginal_power, # sim_corr = sim_corr, # sim_n = n_sim, # test_types = test_types, # test_groups = test_groups, # test_corr = test_corr, # verbose = TRUE # ) # p_sim <- results_calculate_power$details$p_sim # adjusted_p_results <- p_sim # adjusted_significance_level_results <- p_sim # for (j in 1:n_sim) { # temp <- graph_test_closure( # graph, # p_sim[j, ], # alpha = alpha, # test_types = test_types, # test_groups = test_groups, # test_corr = test_corr, # test_values = TRUE # ) # adjusted_p_results[j, ] <- temp$outputs$rejected # intersection <- unique(temp$test_values$results$Intersection) # intersection_rej <- intersection # for (k in 1:length(intersection)) { # intersection_rej[k] <- any(subset( # temp$test_values$results, # Intersection == intersection[k] # )$Inequality_holds) # } # for (k in 1:m) { # id <- substr(intersection, start = k, stop = k) == "1" # adjusted_significance_level_results[j, k] <- suppressWarnings(all(intersection_rej[id])) # } # } # results_adjusted_p <- colMeans(adjusted_p_results) # names(results_adjusted_p) <- names(temp$outputs$rejected) # # Matching power using the adjusted p-value approach # out_adjusted_p[i] <- all.equal( # results_calculate_power$power$power_local, # results_adjusted_p, # tolerance = 1 / n_sim # ) # # results_significance_level <- colMeans(adjusted_significance_level_results) # names(results_significance_level) <- names(temp$outputs$rejected) # # Matching power using the adjusted significance level approach # out_adjusted_significance_level[i] <- all.equal( # results_calculate_power$power$power_local, # results_significance_level, # tolerance = 1 / n_sim # ) # } # out <- cbind(out_adjusted_p, out_adjusted_significance_level) # colnames(out) <- c("adjusted_p", "adjusted_significance_level") # write.csv( # out, # here::here("vignettes/internal-validation_parametric-mixed.csv"), # row.names = FALSE # )