Skip to content

Instantly share code, notes, and snippets.

@gjkerns
Created January 13, 2012 19:30
Show Gist options
  • Save gjkerns/1608265 to your computer and use it in GitHub Desktop.
Save gjkerns/1608265 to your computer and use it in GitHub Desktop.
Simulation for determining sample size in Repeated Measures ANOVA
# 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()))
@ToryMadden
Copy link

Hi gjkerns, I first found this ost of yours on R-bloggers and I'm trying to use it, but I can't run the 'stdevs' command in line 25. Is it supposed to be a standard deviation command or is there something special about it - and if the latter, which package does it belong to? I'm an R and git newbie so I hope you don't mind my asking this question here! -Tory.

@PyRPy
Copy link

PyRPy commented Nov 24, 2020

Hi gjkerns, I first found this ost of yours on R-bloggers and I'm trying to use it, but I can't run the 'stdevs' command in line 25. Is it supposed to be a standard deviation command or is there something special about it - and if the latter, which package does it belong to? I'm an R and git newbie so I hope you don't mind my asking this question here! -Tory.

Hi, Tory
This is a very good example.
'stdevs' is a variable defined in line 9 if you run the code from the beginning. Hope this helps. Good Luck.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment