# R. Chris Fraley | Oct 15, 2003 | http://www.uic.edu/~fraley

# Estimate statistical power for a one way ANOVA

# Plug in hypothesized values, and initialize some variables

# Everything in this section can be set by user

# enter a vector of cell means (one mean per cell | separate means by commas)

means.vec<-c(3.28, 3.07, 4.01)

# enter a vector of cell standard deviations

sd.vec<-c(.90,.82,.71)

# enter the number of people per cell (equal n assumption for simplicity)

n<- 23

# enter the number of simulation trials you would like to run to

# estimate power | test with 50 or so then run full simulation with > 1000

trials<- 100

#-------------------------------------------------

# everything below here can be ignored by the user

#-------------------------------------------------

m<-length(means.vec) # number of cells

# Initialize variables

p.vec.1<-1:trials

xMat<-matrix(0,n,m)

y<-rep(0,n*m)

# Begin loop

for(i in 1:trials){

# Generate data from population values

for(j in 1:m){

xMat[,j]<-rnorm(n,mean=means.vec[j],sd=sd.vec[j])

}

# Find cell means based on simulated sample data

# create vector of dependent values

y<-0

for(j in 1:m){

y<-c(y, xMat[,j])

}

y<-y[-1]

# create fector that codes for iv

a<-0

for(j in 1:m){

a<-c(a, rep(j,n))

}

a<-a[-1]

# Compute p-values based on each F test

results<-summary(aov(y~a))

p.vec.1[i] <-unlist(results)[9]

}

# Print the proportion of trials for which p <= .05

cat("The proportion of trials for which p <= .05","\n")

print(table(p.vec.1<=.05)/trials)