# R/ANOVA_compromise.R In Superpower: Simulation-Based Power Analysis for Factorial Designs

#### Documented in ANOVA_compromise

#' Justify your alpha level by minimizing or balancing Type 1 and Type 2 error rates for ANOVAs.
#' @param design_result Output from the ANOVA_design function
#' @param correction Set a correction of violations of sphericity. This can be set to "none", "GG" Greenhouse-Geisser, and "HF" Huynh-Feldt
#' @param emm Set to FALSE to not perform analysis of estimated marginal means
#' @param emm_model Set model type ("multivariate", or "univariate") for estimated marginal means
#' @param contrast_type Select the type of comparison for the estimated marginal means. Default is pairwise. See ?emmeans::contrast-methods for more details on acceptable methods.
#' @param emm_comp Set the comparisons for estimated marginal means comparisons. This is a factor name (a), combination of factor names (a+b), or for simple effects a | sign is needed (a|b)
#' @param costT1T2 Relative cost of Type 1 errors vs. Type 2 errors.
#' @param priorH1H0 How much more likely a-priori is H1 than H0? Default is 1: equally likely.
#' @param error Either "minimal" to minimize error rates, or "balance" to balance error rates.
#' @param liberal_lambda Logical indicator of whether to use the liberal (cohen_f^2\*(num_df+den_df)) or conservative (cohen_f^2\*den_df) calculation of the noncentrality (lambda) parameter estimate. Default is FALSE.
#' @return Returns dataframe with simulation data (power and effect sizes!), optimal alpha level, obtained beta error rate (1-power/100), and objective (see below for details). If NA is obtained in a alpha/beta/objective columns this indicates there is no effect for this particular comparison. Also returns alpha-beta compromise plots for all comparisons. Note: Cohen's f = sqrt(pes/1-pes) and the noncentrality parameter is = f^2*df(error)
#'
#' \describe{
#'   \item{\code{"aov_comp"}}{A dataframe of ANOVA-level results.}
#'   \item{\code{"aov_plotlist"}}{List of plots for ANOVA-level effects}
#'   \item{\code{"manova_comp"}}{A dataframe of MANOVA-level results.}
#'   \item{\code{"manova_plotlist"}}{List of plots for MANOVA-level effects.}
#'   \item{\code{"emmeans_comp"}}{A dataframe of ANOVA-level results.}
#'   \item{\code{"emm_plotlist"}}{List of plots for estimated marginal means contrasts.}
#'
#' }
#' alpha = alpha or Type 1 error that minimizes or balances combined error rates
#' beta = beta or Type 2 error that minimizes or balances combined error rates
#' objective = value that is the result of the minimization, either 0 (for balance) or the combined weighted error rates
#'
#' @examples
#' \dontrun{
#' design_result <- ANOVA_design(design = "3b*2w",
#' n = 6,
#' mu = c(1, 2, 2, 3, 3, 4),
#' sd = 3,
#' plot = FALSE)
#' example = ANOVA_compromise(design_result,emm = TRUE,emm_comp = "a")
#' }
#' @section References:
#' @importFrom stats optimize
#' @import emmeans
#' @import ggplot2
#' @export
#'

ANOVA_compromise <- function(design_result,
correction = Superpower_options("correction"),
emm = Superpower_options("emm"),
emm_model = Superpower_options("emm_model"),
contrast_type = Superpower_options("contrast_type"),
emm_comp,
costT1T2 = 1,
priorH1H0 = 1,
error = "minimal",
liberal_lambda = Superpower_options("liberal_lambda")){

if (missing(emm_comp)) {
emm_comp = as.character(design_result$frml2)[2] } x1 = ANOVA_exact2(design_result, correction = correction, emm = emm, emm_model = emm_model, emm_comp = emm_comp, verbose = FALSE) aov_comp = data.frame(effect = rownames(x1$main_results),
cohen_f = x1$main_results$cohen_f,
num_df = x1$anova_table$num_df,
den_df = x1$anova_table$den_df,
alpha = NA,
beta = NA,
objective = NA)

if(any(is.na(aov_comp$cohen_f)) || any(is.na(aov_comp$num_df)) ||
any(is.na(aov_comp$den_df)) || any(is.na(aov_comp$effect))) {
stop("Missing aov_comp")
}

aov_plotlist = list()

for (i in 1:nrow(x1$main_results)) { if(x1$main_results[i,]$partial_eta_squared > 1e-11){ run_func = "power.ftest(num_df = x_num_df,den_df = x_den_df,cohen_f = x_cohen_f,alpha_level = x,liberal_lambda = x_lamb)$power/100"
run_func2 = gsub("x_num_df",
aov_comp$num_df[i], run_func) run_func2 = gsub("x_den_df", aov_comp$den_df[i],
run_func2)
run_func2 = gsub("x_cohen_f",
aov_comp$cohen_f[i], run_func2) run_func2 = gsub("x_lamb", liberal_lambda, run_func2) alpha_res = optimal_alpha(power_function = run_func2, costT1T2 = costT1T2, priorH1H0 = priorH1H0, error = error, plot = FALSE) aov_comp[i,]$alpha = alpha_res$alpha aov_comp[i,]$beta = alpha_res$beta aov_comp[i,]$objective = alpha_res$objective aov_name = aov_comp[i,]$effect
alpha_res$plot = alpha_res$plot + ggtitle(aov_name)
aov_plotlist[[aov_name]] = alpha_res$plot } } if (!is.null(x1$manova_results)) {
manova_comp = data.frame(effect = rownames(x1$manova_results), cohen_f = x1$manova_results$cohen_f, num_df = x1$manova_table$num_df, den_df = x1$manova_table$den_df, alpha = NA, beta = NA, objective = NA) manova_plotlist = list() for (i in 1:nrow(manova_comp)) { if (x1$manova_table[i,]$p.value < 1) { run_func = "power.ftest(num_df = x_num_df,den_df = x_den_df,cohen_f = x_cohen_f,alpha_level = x,liberal_lambda=x_lamb)$power/100"
run_func2 = gsub("x_num_df",
manova_comp$num_df[i], run_func) run_func2 = gsub("x_den_df", manova_comp$den_df[i],
run_func2)
run_func2 = gsub("x_cohen_f",
manova_comp$cohen_f[i], run_func2) run_func2 = gsub("x_lamb", liberal_lambda, run_func2) alpha_res = optimal_alpha( power_function = run_func2, costT1T2 = costT1T2, priorH1H0 = priorH1H0, error = error, plot = FALSE ) manova_comp[i,]$alpha = alpha_res$alpha manova_comp[i,]$beta = alpha_res$beta manova_comp[i,]$objective = alpha_res$objective manova_name = manova_comp[i,]$effect
alpha_res$plot = alpha_res$plot + ggtitle(manova_name)
manova_plotlist[[manova_name]] = alpha_res$plot } } manova_comp = data.frame(effect = manova_comp$effect,
cohen_f = manova_comp$cohen_f, num_df = manova_comp$num_df,
den_df = manova_comp$den_df, alpha = manova_comp$alpha,
beta = manova_comp$beta, objective = manova_comp$objective)

} else {
manova_comp = NULL
manova_plotlist = NULL
}

if (!is.null(x1$emm_results)) { emmeans_comp = data.frame(effect = x1$emm_results$contrast, cohen_f = x1$emm_results$cohen_f, num_df = 1, den_df = x1$emmeans_table$df, alpha = NA, beta = NA, objective = NA) emm_plotlist = list() for (i in 1:nrow(x1$emm_results)) {
if (x1$emmeans_table[i,]$p.value <  1) {
run_func = "power.ftest(num_df = x_num_df,den_df = x_den_df,cohen_f = x_cohen_f,alpha_level = x,liberal_lambda=x_lamb)$power/100" run_func2 = gsub("x_num_df", emmeans_comp$num_df[i],
run_func)
run_func2 = gsub("x_den_df",
emmeans_comp$den_df[i], run_func2) run_func2 = gsub("x_cohen_f", emmeans_comp$cohen_f[i],
run_func2)
run_func2 = gsub("x_lamb",
liberal_lambda,
run_func2)
alpha_res = optimal_alpha(
power_function = run_func2,
costT1T2 = costT1T2,
priorH1H0 = priorH1H0,
error = error,
plot = FALSE
)
emmeans_comp[i,]$alpha = alpha_res$alpha
emmeans_comp[i,]$beta = alpha_res$beta
emmeans_comp[i,]$objective = alpha_res$objective
emm_name = emmeans_comp[i,]$effect alpha_res$plot = alpha_res$plot + ggtitle(emm_name) emm_plotlist[[emm_name]] = alpha_res$plot
}
}
emmeans_comp = data.frame(effect = emmeans_comp$effect, cohen_f = emmeans_comp$cohen_f,
num_df = emmeans_comp$num_df, den_df = emmeans_comp$den_df,
alpha = emmeans_comp$alpha, beta = emmeans_comp$beta,
objective = emmeans_comp\$objective)

} else {
emmeans_comp = NULL
emm_plotlist = NULL
}

# S3 method
#class(compromise_ANOVA) <- "compromise_ANOVA"
#attr(compromise_ANOVA, aov_comp)

structure(list(aov_comp = aov_comp,
aov_plotlist = aov_plotlist,
manova_comp = manova_comp,
manova_plotlist = manova_plotlist,
emmeans_comp = emmeans_comp,
emm_plotlist = emm_plotlist,
method = "ANOVA_compromise"
),
class = "opt_alpha")

}


## Try the Superpower package in your browser

Any scripts or data that you put into this service are public.

Superpower documentation built on May 17, 2022, 5:08 p.m.