# R. Chris Fraley | Oct 15, 2001 | updated Nov 29 2002
# Estimate statistical power for a between-subjects 2 x 2 ANOVA
 
# To use this program, simply plug in the expected or hypothesized means
# for each of the four cells in the 2 by 2 design.  
# The means are represented as m1, m2, m3, and m4 below.
# By default, I've set these to 10.  
#
# Step 1: Please change the means (e.g., m1, m2) to values that
# are appropriate for your research example.
# Step 2: Set k to equal the hypothesized population standard deviation.
# Step 3: Set n to the number of subjects in each cell
#
# At the end of the simulation, the program will tell you the proportion of times you'll
# obtain a p-value less than 5% for the main effect of x, main effect of y, and the interaction 
# of x and y.
 
 
#                   variable y
# variable x    level 1    level 2
#             -----------------------    
#    level 1 |  cell 1  |   cell 2  |
#            ------------------------ 
#    level 2 |  cell 3  |   cell 4  |
#            -----------------------
 
m1<- 10    # hypothesized mean cell 1
m2<- 10    # hypothesized mean cell 2
m3<- 10    # hypothesized mean cell 3
m4<- 10    # hypothesized mean cell 4
 
k<-2         # hypothesized sd within each condition
n<-10        # number of people per cell
trials<-1000 # number of trials for simulation
 
#-----------------------------------------------------------
# The rest of this doesn't need to be adjusted by the user.
#-----------------------------------------------------------
 
m<-4         # number of cells 
 
p.vec.1<-1:trials   
p.vec.2<-1:trials
p.vec.1x2<-1:trials
 
# Begin loop
 
for(i in 1:trials){
 
 
# Generate data from population values
 
  x1<-rnorm(n,mean=m1,sd=k) # simulate data for cell 1
  x2<-rnorm(n,mean=m2,sd=k) # simulate data for cell 2
  x3<-rnorm(n,mean=m3,sd=k) # simulate data for cell 3
  x4<-rnorm(n,mean=m4,sd=k) # simulate data for cell 4
 
# Find cell means based on simulated sample data
 
y<-c(x1,x2,x3,x4)
a<-c(rep(0,n*2),rep(1,n*2))
b<-c(rep(0,n),rep(1,n),rep(0,n),rep(1,n))
a<-as.factor(a)
b<-as.factor(b)
 
# Compute p-values based on each F test
 
results<-summary(aov(y~a*b))
p.vec.1[i]  <-results[1,5]
p.vec.2[i]  <-results[2,5]
p.vec.1x2[i]<-results[3,5]
 
}
 
# Print the proportion of p<.05 for each effect

cat("Proportion of times the null hypothesis will be rejected for a main effect of the first factor ","\n") 
print(table(p.vec.1<=.05)/trials)
cat("Proportion of times the null hypothesis will be rejected for a main effect of the second factor ","\n") 
print(table(p.vec.2<=.05)/trials)
cat("Proportion of times the null hypothesis will be rejected for the interaction ","\n") 
print(table(p.vec.1x2<=.05)/trials)