Power simulation for mixed design ANOVA

anova.pwr.mixed.sim.r

require(mvtnorm)
# function --------------------------------------------------------------------
anova.pwr.mixed.sim <- function(data, Formula, FactorA, n, rho, sims) {
  Formula <- formula(Formula) # convert to formula
  m <- nlevels(data$time) # number of time points
  k <- nlevels(data[[FactorA]]) # number of groups
  if(length(n) == 1) n <- rep(n, k) # n for each level of 'group'
  n2 <- rep(n, each=m) # n for each cell
  n_index <- cumsum(n)  # help-variable used when creating subjects factor
  subjects <- NULL
  # create subjects factor
  for(i in 1:k) { 
    if(i==1) subjects <- c(subjects, (rep(1:n_index[i], m)))    
    else subjects <- c(subjects, rep((n_index[i-1]+1):n_index[i], m))
  }
  subjects <- factor(subjects)
  # expand model in 'dataframe' to fit subjects
  dataframe <- NULL 
  for(i in 1:nrow(study)) {  
    dataframe <- rbind(dataframe, data[rep(i, n2[i]),])
  }
  dataframe$subjects <- subjects
  
  # create covariance matrix from rho and sigmas within cells
  cov_matrix <- function(sigmas) {
    B <- matrix(sigmas, ncol=length(sigmas), nrow=length(sigmas))
    sigma <- t(B) * B * rho
    diag(sigma) <- sigmas^2 
    sigma
  }
  mu_m <- matrix(data$DV, ncol=m, byrow=T) # used to generate data 
  sigma_m <-  matrix(data$SD, ncol=m, byrow=T) # used to generate data
  nterms <- 3 # number of possible effects
  sig <- matrix(ncol=nterms, nrow=sims) # pre allocate sig matrix
  # actual simulation
  for (i in 1:sims) {
    dataframe$DV <- unlist(lapply(1:k, function(i) rmvnorm(n[i], mu_m[i,], cov_matrix(sigma_m[i,])))) # generate data for all cells
    result <- summary(aov(Formula, dataframe)) # perform ANOVA
    tmp <- c(result[[1]][[1]][[5]], result[[2]][[1]][[5]]) # get p-values
    tmp_matrix <- matrix(tmp[!is.na(tmp)], nrow=1) # remova NA
    sig[i,] <- tmp_matrix # save p-values
  }
  terms_b <- head(rownames(result[[1]][[1]]), 1) # get names of between-subjects terms   
  terms_w <- head(rownames(result[[2]][[1]]), 2)  # get names of within-subjects terms
  out <- data.frame("terms" = c(terms_b, terms_w), "power" =  colMeans(sig < 0.05)) # output
  out # return output
}