public
Last active

Simulation for determining sample size in Repeated Measures ANOVA

  • Download Gist
powerSampleSize.R
R
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60
# Simulation study for sample size between/within
# got Treat + Sham between subjects
# got Time within subjects
 
nPerGroup <- 30
nTime <- 4
muTreat <- c(37, 32, 20, 15)
muSham <- c(37, 32, 25, 22)
stdevs <- c(12, 10, 8, 6)
stdiff <- 9
nSim <- 1000
 
# 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
ones <- rep(1, nTime)
A <- stdevs^2 %o% ones
B <- (A + t(A) + (stdiff^2)*(diag(nTime) - ones %o% ones))/2
 
 
# can run it through once to check that it works
library(MASS)
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)
summary(aovComp)
 
# 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()))

Please sign in to comment on this gist.

Something went wrong with that request. Please try again.