Created

Embed URL

HTTPS clone URL

SSH clone URL

You can clone with HTTPS or SSH.

Download Gist

Simulation for determining sample size in Repeated Measures ANOVA

View powerSampleSize.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()))
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Something went wrong with that request. Please try again.