R/exact_twoway_anova_power.R

In extraSuperpower: Power Calculation for Two-Way Factorial Designs

#' Two-way factorial ANOVA exact sample size calculation for independent samples
#'
#' This functions takes the effect sizes (Cohen's f) for two main effects and their interaction and estimates power
#' a range of sample sizes. The input for this function can be generated by `effsize`.
#'
#' @param a Number of levels of the first factor
#' @param b Number of levels of the second factor
#' @param effect_sizes Numeric vector of length 3. The first two elements are the effect sizes for the main effects of the first
#' and second factors, respectively. The third element is the interaction effect size.
#' @param n Number of experimental units in each group for which power (1-beta) will be calculated.
#' @param alpha numeric - Type I error probability. Defaults to 0.05
#' @param factor_names character - vector of length 2. Names of the 2 factors to be evaluated. Default is to inherit names
#' from effect_sizes. If effect_sizes has no names and no factor_names are provided, factors will be named 'FactorA' and
#' 'FactorB'.
#' @param plot logical - Should the power curve be plotted. Default is TRUE.
#' @param target_power numeric - Desired power to be attained. Accepts values between 0 and 1, defaults to 0.8.
#' @param title character - Title for the graph. Defaults to 'Power curve from exact ANOVA test'
#' @param target_line logical - Set to TRUE. If FALSE no target line will be drawn. Overrides target_power.
#' @param alpha_line logical - Should a line at the set type I error be plotted
#'
#' @details Probably the best way to calculate power for independent balanced designs
#'
#' @return A list that contains the number of levels for each factor, the chosen significance level and a data.frame in which the
#' first column is the group sample size and the remaining three columns are the power for the main effect of the first
#' factor, the main effect of the second factor and their interaction, respectively.
#'
#' @return Optionally, a graph that displays the power curves.
#'
#' @importFrom stats qf
#' @importFrom stats pf
#'
#' @examples
#' refmean <- 1
#' treatgroups <- 4
#' timepoints <- 5
#' treateff <- 1.25
#' timeeff <- 0.85
#' cellswithinteraction <- matrix(c(rep(2,3), 3:5), 3,2)
#' #second level of factor A interacts with 3rd, 4th and 5th level of factor B
#'
#'factors_levels_names <- list(treatment=letters[1:treatgroups], time=1:timepoints)
#'
#' effects_treat_time_interact <- calculate_mean_matrix(refmean = refmean,
#'                                                      nlfA = treatgroups, nlfB = timepoints,
#'                                                      fAeffect = treateff, fBeffect = timeeff,
#'                                                      label_list = factors_levels_names,
#'                                                      groupswinteraction = cellswithinteraction,
#'                                                      interact=1.3)
#'
#' fxs <- effsize(effects_treat_time_interact)
#' exact_twoway_anova_power(a= treatgroups, b=timepoints, effect_sizes=fxs, n=5:20)
#'
#' @export
exact_twoway_anova_power <- function(a, b, effect_sizes, n, alpha = 0.05, factor_names = NULL, plot=TRUE, title = NULL, target_power=NULL, target_line=TRUE, alpha_line=TRUE)
{
  if (length(a)!=1 | length(b)!=1 | a %% 1 != 0 | b %% 1 != 0 | a == 0 | b == 0)
  {
    stop("Number of factor levels a and b must be single positive integers")
  }
  if (length(effect_sizes)!=3 | any(effect_sizes<0) | !is.numeric(effect_sizes))
  {
    stop("Effect sizes must be a vector of 3 numbers with value greater or equal to zero")
  }
  if (any(n %% 1 != 0) | any(n==0))
  {
    stop("n must be a vector of positive integers (may be a vector of length 1)")
  }
  if (!is.null(factor_names) & length(factor_names) != 2)
  {
    stop("factor_names must be a vector of length 2")
  }
  calcpow <- function(num_df, f_size)
  {
    lambda <- N * f_size^2
    q <- qf(alpha, num_df, denom_df, lower.tail = FALSE)
    power <- pf(q, num_df, denom_df, lambda, lower.tail = FALSE)
    return(power)
  }
  fAdf <- a - 1
  fBdf <- b - 1
  ABdf <- (a-1)*(b-1)
  powercurve <- NULL
  ##go through possible sample sizes, calculate power
  for (i in n)
  {
    N <- a*b*i
    denom_df <- a*b*(i-1)
    fApow <- calcpow(fAdf, effect_sizes[1])
    fBpow <- calcpow(fBdf, effect_sizes[2])
    inTpow <- calcpow(ABdf, effect_sizes[3])
    res <- cbind(i, fApow, fBpow, inTpow)
    powercurve <- rbind(powercurve, res)
  }
  if(!is.null(factor_names))
  {
    factor_names <- c(factor_names, paste(factor_names, collapse = ":"))
  } else if(!is.null(colnames(effect_sizes)) & length(colnames(effect_sizes))==3)
  {
    factor_names <- colnames(effect_sizes)
  } else
  {
    factor_names <- c("FactorA", "FactorB", "Interaction")
  }
  powercurve <- data.frame(powercurve)
  names(powercurve) <- c("n", factor_names)
  longpowercurve <- reshape2::melt(powercurve, id.var="n", value.name = "power", variable.name="effect")
  NOTE <- paste("n is number in each group, total sample = n *",
                a * b)
  METHOD <- "Balanced two-way analysis of variance power calculation"
  if(!plot)
  {
    return(list(a = a, b = b, sig.level = alpha, powercurve = powercurve, note=NOTE,
       method = METHOD))
  } else if(plot)
    exactpower <- list(a = a, b = b, sig.level = alpha, powercurve = powercurve, note=NOTE,
                       method = METHOD)
   list(exactpower = exactpower, powerplot = plot_powercurves(longpowercurve))
}