# 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)