Skip to content

Instantly share code, notes, and snippets.

Embed
What would you like to do?
Analytical power analysis for mixed design ANOVA
# 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
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment