#This code has been modified from Aguilar-Sanchez (2009) and Gelman and Hill (2007)
#Install the "arm" package from the CRAN

#load the "arm" package

library("arm")

#This function createss a set of variables representative of my data (the entire data set)

estar.fake <- function (J, K){
  speaker <- rep(1:J, each=K)       # creates a vector of numbers from 1 to J to represent each speaker
  subject <- sample(c(1:0), K, replace=T, prob=NULL) #creates a vector of 0s and 1s representing the variable subject
  adjectiveclass <- sample(c(1:0), K, replace=T, prob=NULL) #creates a vector of categories from 1 to 6 representing the variable Adjective Class
  understructure <- sample(c(1:5), K, replace=T, prob=NULL) #creates a vector of categories from 1 to 5 respresenting the variable underlying structure
  gradiency <- sample(c(1:0), K, replace=T, prob=NULL)    #createes a vector of 0s and 1s representing the variable gradiency 
  predicatetype <- sample(c(1:0), K, replace=T, prob=NULL) #createes a vector of 0s and 1s representing the variable predicate reading 
  succhange <- sample(c(1:0), K, replace=T, prob=NULL)   #createes a vector of 0s and 1s representing the variable succeptibility to change 
  frameref <- sample(c(1:0), K, replace=T, prob=NULL)  #createes a vector of 0s and 1s representing the variable frame of reference 
  expreferent <- sample(c(1:3), K, replace=T, prob=NULL)  #createes a vector of categories from 1 to 3 representing the variable experience with the referent.
  resultstate <- sample(c(1:0), K, replace=T, prob=NULL) #createes a vector of 0s and 1s representing the variable resultant state 
  age <- rnorm(J, 40, 5)		#creates a vector of ages with mean of 40 and SD of 5
  age1 <- age[speaker]			#converts the age vector into a predictor variable
  gender <- sample(c(1:0), J, replace=T, prob=NULL)     #creates a vector of 0s and 1s representing the variable gender
  gender1 <- gender[speaker]		#creates a vector 0s and 1s of gender for each column according to the speaker
  educ <- rnorm (J, 12, 3)		#creates a vector of normally distributed data with mean of 12 and SD of 3 representing years of education
  educ1 <- educ[speaker]		#creates the predictor variable education
  biling <- sample (c(1:0), J, replace=T, prob=NULL)      #creates a vector of 0s and 1s representing the variable bilinguism with a probability of x=1=.8
  biling1 <- biling[speaker]		#creates a predictor variable called bilinguism
#				    #now we create the vector of 1s and 0s representing the Dependent variable
estar <- sample(c(1:0), K, replace=T, prob=NULL)
return (data.frame (estar, subject, adjectiveclass, understructure, gradiency, predicatetype, succhange, frameref, expreferent, resultstate, speaker, gender1, age1, educ1, biling1))
}

save.image("/Users/user/Desktop/simulation_10_04_11")  #erase this or modify it to match your computer's needs

#This function to test for power for the variable coefficient of the model without gradiancy, underlinestructure, and bilingualism variables. 
#The number 1000 in this function is essential to run the Monte Carlo Simulations

estar.power <- function (J, K, n.sims=1000) {
  signif <- rep (NA, n.sims)
  for (s in 1:n.sims) {
    estar.data <- estar.fake (J, K)   
 lme.power <- lmer(estar ~ subject + adjectiveclass + predicatetype + succhange + frameref + expreferent + resultstate + (1 |gender1) + (1| age1) + (1| educ1) , family=binomial(link="logit"), data=estar.data) #I eliminated the variable bilingualism from the model to exemplify the effects of using less social variables
  theta.hat <- fixef(lme.power)
  theta.se <- se.fixef(lme.power)
  signif[s] <- ifelse (sqrt(theta.hat^2) - (2*theta.se) > 0,1,0)
}
power <- mean (signif)
return (power)
}

estar.power(J=60, K=30, n.sims=10)  # this code tests the function estar.power

save.image("/Users/user/Desktop/simulation_10_04_11")  #erase this or modify it to match your computer's needs

#Monte Carlo simulation
#Loop to get the power calculations

J.values<- c(20, 30, 60, 100, 200, 300)
n.sims.values <- rep(100,9)    #change this value back to 1000 in order for it to be correct
K.values <- c(25, 30, 35)
power.values <- array (NA, c(length(J.values),length(K.values)))
for (i1 in 1:length(J.values)){
  for (i2 in 1:length(K.values)){
    cat ("computing power calculation for J =", J.values[i1], ", K =", K.values[i2], "\n")
    power.values[i1,i2] <- estar.power (J=J.values[i1], K=K.values[i2], n.sims=n.sims.values[i1])
    cat ("power =", power.values[i1,i2], "\n")
  }
}


#This code will give you the plots needed to determine the sample size you need for your study
#plot all the cures

plot (c(0,max(J.values)), c(0,1), xaxs="i", yaxs="i", xlab="number of speakers", ylab="power", type="n")
  for (i2 in 1:length(K.values)){
    lines (c(0,J.values), c(.025, power.values[,i2]))
}

save.image("/Users/user/Desktop/simulation_10_07_11")

#10-05-10 Started at 8:45 a.m.  and Ended at 8:57 a.m. 10-06-11  #erase this or modify it to match your computer's needs (windows use "\" and MACs use "/"