R version 2.13.2 (2011-09-30)
Copyright (C) 2011 The R Foundation for Statistical Computing
ISBN 3-900051-07-0
Platform: i386-pc-mingw32/i386 (32-bit)
R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.
Natural language support but running in an English locale
R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.
Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.
> local({pkg <- select.list(sort(.packages(all.available = TRUE)),graphics=TRUE)
+ if(nchar(pkg)) library(pkg, character.only=TRUE)})
Loading required package: MASS
Loading required package: Matrix
Loading required package: lattice
Attaching package: 'Matrix'
The following object(s) are masked from 'package:base':
det
Loading required package: lme4
Attaching package: 'lme4'
The following object(s) are masked from 'package:stats':
AIC, BIC
Loading required package: R2WinBUGS
Loading required package: coda
Attaching package: 'coda'
The following object(s) are masked from 'package:lme4':
HPDinterval
Loading required package: abind
arm (Version 1.4-13, built: 2011-6-19)
Working directory is C:/Users/aguilar.jorg/Documents
Loading required package: foreign
Attaching package: 'arm'
The following object(s) are masked from 'package:coda':
traceplot
> 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 (windows use "\" and MACs use "/"
>
> #This function to test for power of the saturated model
> #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 + understructure + gradiency + predicatetype + succhange + frameref + expreferent + resultstate + (1 |gender1) + (1| age1) + (1| educ1) + (1| biling1), family=binomial(link="logit"), data=estar.data)
+ theta.hat <- fixef(lme.power) #this represents the coefficient Beta yielded by the lmer() function
+ theta.se <- se.fixef(lme.power) #this represents the standard error for the saturated model yielded by the lmer() function
+ signif[s] <- ifelse (sqrt(theta.hat^2) - (2*theta.se) > 0,1,0) #this code tests the significance at the .05 level or 2 se above or below the mean
+ }
+ power <- mean (signif)
+ return (power)
+ }
>
> estar.power(J=60, K=30, n.sims=10) # this code tests the function estar.power and only runs the estar.power function 10 times
[1] 0.8
There were 13 warnings (use warnings() to see them)
>
> 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(10,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")
+ }
+ }
computing power calculation for J = 20 , K = 25
power = 0.7
computing power calculation for J = 20 , K = 30
power = 0.8
computing power calculation for J = 20 , K = 35
power = 0.8
computing power calculation for J = 30 , K = 25
power = 0.3
computing power calculation for J = 30 , K = 30
power = 0.7
computing power calculation for J = 30 , K = 35
power = 0.8
computing power calculation for J = 60 , K = 25
power = 0.6
computing power calculation for J = 60 , K = 30
power = 0.9
computing power calculation for J = 60 , K = 35
power = 0.7
computing power calculation for J = 100 , K = 25
power = 0.4
computing power calculation for J = 100 , K = 30
power = 0.8
computing power calculation for J = 100 , K = 35
power = 0.9
computing power calculation for J = 200 , K = 25
power = 0.8
computing power calculation for J = 200 , K = 30
power = 1
computing power calculation for J = 200 , K = 35
power = 0.9
computing power calculation for J = 300 , K = 25
power = 0.9
computing power calculation for J = 300 , K = 30
power = 0.8
computing power calculation for J = 300 , K = 35
power = 1
There were 50 or more warnings (use warnings() to see the first 50)
>
>
> #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]))
+ }
>