ex power simulation for repeated measures

ex power sim repeated measures anova.R

# *-----------------------------------------------------------------
# | PROGRAM NAME: ex power sim repeated measures anova.R 
# | DATE: 10/10/20 
# | CREATED BY: MATT BOGARD 
# | PROJECT FILE: https://gist.github.com/BioSciEconomist           
# *----------------------------------------------------------------
# | PURPOSE: annotated code based on reference cited below
# *----------------------------------------------------------------

# reference: https://www.r-bloggers.com/2012/01/power-and-sample-size-for-repeated-measures-anova-with-r/
# Github: https://gist.github.com/gjkerns/1608265  

# Simulation study for sample size between/within
# got Treat + Sham between subjects
# got Time within subjects

nPerGroup <- 30 # sample size per group
nTime     <- 4  # number of time periods or measures
muTreat   <- c(37, 32, 20, 15) # average outcome values for the treatment
muSham    <- c(37, 32, 25, 22) # average outcome values for control
stdevs    <- c(12, 10, 8, 6)   # standard deviations 
# assumptions - from other studies indicate NDI scores have a standard deviation
# of around 12, and those have been observed to decrease over time.
stdiff    <- 9  # standard deviation of differences sphericity assumes constant          #
nSim      <- 1000 # number of simulations

# set up the indep var data
Subject <- factor(1:(nPerGroup*2))
Time <- factor(1:nTime, labels = c("0min", "15min", "48hrs", "96hrs"))

theData <- expand.grid(Time, Subject)
names(theData) <- c("Time", "Subject")

tmp <- rep(c("Treat", "Sham"), each = nPerGroup * nTime)
theData$Method <- factor(tmp)

# to set up variance-covariance matrix
# based on assumptions cited in:
# Huynh Huynh, and Leonard S. Feldt. “Conditions Under Which Mean Square Ratios
# in Repeated Measurements Designs Have Exact F-Distributions.” Journal of 
# the American Statistical Association, vol. 65, no. 332, 1970, pp. 1582–1589. 
# JSTOR, www.jstor.org/stable/2284340. Accessed 8 Oct. 2020.

ones <- rep(1, nTime) # vector of ones 
A <- stdevs^2 %o% ones 
B <- (A + t(A) + (stdiff^2)*(diag(nTime) - ones %o% ones))/2

# note - the %o% is the outer product operator 
1:3 %o% 1:3 # example

# can run it through once to check that it works
library(MASS) # library used to simulate multivariate normal distributions
tmp1 <- mvrnorm(nPerGroup, mu = muTreat, Sigma = B) # simulate treatment outcome
tmp2 <- mvrnorm(nPerGroup, mu = muSham, Sigma = B) # simulate control outcome
theData$NDI <- c(as.vector(t(tmp1)), as.vector(t(tmp2))) # merge into dataframe
aovComp <- aov(NDI ~ Time*Method + Error(Subject/Time), theData) # fit model run
summary(aovComp) # results

# some descriptive statistics and graphs
print(model.tables(aovComp, "means"), digits = 3) 
boxplot(NDI ~ Time, data = theData)              
boxplot(NDI ~ Method, data = theData)            
boxplot(NDI ~ Time*Method, data = theData)       
with(theData, interaction.plot(Time, Method, NDI))

###############################################
# for power estimate run the below
# don't forget to set up theData and var-cov
library(MASS)

runTest <- function(){
  tmp1 <- mvrnorm(nPerGroup, mu = muTreat, Sigma = B)
  tmp2 <- mvrnorm(nPerGroup, mu = muSham, Sigma = B)
  theData$NDI <- c(as.vector(t(tmp1)), as.vector(t(tmp2)))
  aovComp <- aov(NDI ~ Time*Method + Error(Subject/Time), theData) 
  b <- summary(aovComp)$'Error: Subject:Time'[[1]][2,5]
  b < 0.05
}

# here is estimate of power for given nPerGroup
mean(replicate(nSim, runTest()))