Analytical power analysis for mixed design ANOVA

anova.pwr.mixed.r

# function ---------------------------------------------------------------
anova.pwr.mixed <- function(data, Formula, n, m, rho, n.needed=FALSE) {
  require(stringr)
  Formula <- formula(paste(Formula, "+Error(group)"))               # aov formula
  x <- summary(aov(Formula, data))          # aov model
  terms <- c(rownames(x[[1]][[1]]), rownames(x[[2]][[1]])) # terms in model
  terms <- str_pad(terms, max(nchar(terms)), "right")      # fix typography
  
  pwr <- function(n) {
    k <- length(unique(study$group))          # n groups
    sigma <- mean(data$SD)                    # sigma
    data$n <- n                               # n per group
    N <- sum(with(data, tapply(n,group,unique))) # total N
    x <- summary(aov(Formula, data))          # aov model
    var_b <- x[[1]][[1]][[2]] / (k*m)         # SS between 
    var_w <- x[[2]][[1]][[2]] / (k*m)         # SS within
    f_between <- sqrt(var_b / sigma^2)        # f between
    f_within <- sqrt(var_w / sigma^2)         # f within
    
    # Degrees of freedom ------------------------------------------------------
    df1_b <- x[[1]][[1]][[1]]              # df model between
    df1_w <- x[[2]][[1]][[1]]              # df model within
    df1 <- c(df1_b, df1_w)
    df2_b <- N-1-sum(df1_b)                # df error between
    df2_w <- (N*m-1) - (N-1) - sum(df1_w)  # df error within
    
    # Power Between -----------------------------------------------------------------
    lambda <- N * (m/(1+(m-1)*rho)) * f_between^2        # lambda between
    F_critical <- qf(0.05,df1_b,df2_b, lower.tail=FALSE) # F critical between
    power_b <- pf(F_critical, df1_b, df2_b, lambda, lower.tail=FALSE) # power between
    
    # Power Within -----------------------------------------------------------------
    lambda <- (N * m * f_within^2) / (1 - rho)           # lambda within
    F_critical <- qf(0.05,df1_w,df2_w, lower.tail=FALSE) # F critical within
    power_w <- pf(F_critical, df1_w, df2_w, lambda, lower.tail=FALSE) # power within
    c("b" = power_b, "w" = power_w)
  }
  
  # Get n for desired power -------------------------------------------------
  power <- pwr(n) # get powr for 'n'
  if(n.needed) {
    n_needed <- rep(NA, length(power)) # pre-create object
    # loop til desired power is found
    for(i in 1:length(power)) {
      n2 <- 5
      while(pwr(n2)[i] < 0.8) n2 <- n2+1
      n_needed[i] <- n2 
    }
  } else n_needed <- NA

  
  # output ------------------------------------------------------------------
  out <- data.frame(terms, "power" = round(power, 3),
                    "n.needed" = n_needed) # n per group
  names(out)[1] <- str_pad("Terms", max(nchar(terms)), "right")
  out
}