# Power analyses (11-01-2008) for Boundary N+1/N+2

# use estimates from the lmer analyses (log transformed data):
       # (Intercept)           x1            x2         x3     freqclow
# FFD n      5.388685  0.003336719 -0.0142235653 0.01170662 7.739181e-05
# FFD n+1    5.488482 -0.011530990 -0.0001957424 0.15032796 6.515374e-02
# FFD n+2    5.422299 -0.009868352  0.0186719449 0.04606185 4.188075e-02
# GZD n      5.552655  0.016749127 -0.0251055874 0.04432217 1.000379e-02
# GZD n+1    5.613172 -0.004258721  0.0276584508 0.25704613 1.262018e-01
# GZD n+2    5.536700 -0.008882073  0.0118849734 0.03326725 6.596512e-02
          # s.sn      s.id       s.y    N
# FFD n   0.06394302 0.1129541 0.2574435 4228
# FFD n+1 0.04461764 0.1382883 0.2953375 4296
# FFD n+2 0.07860019 0.1361069 0.2824953 4239
# GZD n   0.13545470 0.1805066 0.3467311 4227
# GZD n+1 0.07148875 0.1669347 0.3290196 4289
# GZD n+2 0.12981321 0.1558093 0.3403959 4238
# Power analyses using simulated data (Gelman & Hill, 2007)
#
rm(list=ls())

workingdir <- choose.dir()

#Alternative: Set working directory directly:

# workingdir <- "/working directory"
setwd(workingdir)

# directory must contain "lmer estimates for power.log.rData" (generated by "lmer Analyses.R")
# as well as the script files "ffd.fake.log.R" and "n2power.R"

load(file="lmer estimates for power analysis log.rData") # contains all.f

require(lme4)
require(reshape)
source("ffd.fake.log.R")

source("n2power.R")

# Angele et al. (K=152, J=32)

K <- 152
J <-  32

########N.SIMS#############
# Large numbers of simulations (1000 and up) are desirable, but can take several hours to complete on an average machine.
# Test the script with few simulations (e.g. 10), then increase the number of simulations once everything works.

n.sims <- 10



##########################

#write results to file
sink("N12power.log.txt")
# simulations for FFD
for(i in 1:3)
{

# use real Ns from the analyses to estimate proportion of missing data
prb1 <- all.f[i,"N"]/4712

cat("\n")
# first do a simulation assuming no data loss
cat(paste("\nProp = 100, n.sims =", n.sims,"\n"))

cat(paste(rownames(all.f)[i],"\n"))
print(all.f[i,])
cat("\n\n")

ptm <- proc.time()
# assuming that n+2 effects have an effect size of 7 ms
ffd.pow100 <- n2.power(J, K, prop=1,   n.sims=n.sims, effect.size=7, mu.a.true = 200, sigma.item.true = all.f[i,"s.sn"], sigma.id.true = all.f[i,"s.id"], sigma.y.true = all.f[i,"s.y"] )     #  effect.size = 7
# how long do the simulations take to complete?
print(proc.time() - ptm)

cat("\n")
#print probability of detecting an effect of 7 ms, i.e. power

print(cast(melt(as.data.frame(ffd.pow100)), variable~., mean))
#now do the same simulations using the correct proportion of missing data
cat(paste("\nProp =", round(prb1, digits=2)," n.sims =", n.sims,"\n"))

cat("\n")
cat(paste(rownames(all.f)[i],"\n"))
print(all.f[i,])
cat("\n\n")


ptm <- proc.time()

ffd.pow60 <- n2.power(J, K, prop=prb1,   n.sims=n.sims, effect.size=7, mu.a.true = 230, sigma.item.true = all.f[i,"s.sn"], sigma.id.true = all.f[i,"s.id"], sigma.y.true = all.f[i,"s.y"] )     #  effect.size = 7

print(proc.time() - ptm)

cat("\n")


#print probability of detecting an effect of 7 ms, i.e. power
print(cast(melt(as.data.frame(ffd.pow60)), variable~., mean))
save(list = c("ffd.pow100","ffd.pow60"), file = paste("rpow_raw_",i,".RData",sep=""))
}

# now do the same for GD
for(i in 4:6)
{

prb1 <- all.f[i,"N"]/4712

cat("\n")

cat(paste("\nProp = 100, n.sims =", n.sims,"\n"))

cat(paste(rownames(all.f)[i],"\n"))
print(all.f[i,])
cat("\n\n")

ptm <- proc.time()

ffd.pow100 <- n2.power(J, K, prop=1,   n.sims=n.sims, effect.size=7, mu.a.true = 200, sigma.item.true = all.f[i,"s.sn"], sigma.id.true = all.f[i,"s.id"], sigma.y.true = all.f[i,"s.y"] )     #  effect.size = 7

print(proc.time() - ptm)

cat("\n")

#for (i in 1:4) print(mean(ffd.pow100[,i]))            #  .977,  .979,   .476

print(cast(melt(as.data.frame(ffd.pow100)), variable~., mean))

cat(paste("\nProp =", round(prb1, digits=2)," n.sims =", n.sims,"\n"))

cat("\n")
cat(paste(rownames(all.f)[i],"\n"))
print(all.f[i,])
cat("\n\n")


ptm <- proc.time()

ffd.pow60 <- n2.power(J, K, prop=prb1,   n.sims=n.sims, effect.size=7, mu.a.true = 200, sigma.item.true = all.f[i,"s.sn"], sigma.id.true = all.f[i,"s.id"], sigma.y.true = all.f[i,"s.y"] )     #  effect.size = 7

print(proc.time() - ptm)

cat("\n")

#for (i in 1:4) print(mean(ffd.pow60[,i]))             #  .865,  .871,   .858

print(cast(melt(as.data.frame(ffd.pow60)), variable~., mean))
save(list = c("ffd.pow100","ffd.pow60"), file = paste("rpow_raw_log_",i,".RData",sep=""))
}


sink(file=NULL)

# Power estimates from the simulations are in N12power.log.txt
# e.g:
   # variable (all)
# 1  signif.f   0.7 Frequency effect
# 2 signif.x1   0.4 Effect of N+2 preview when N+1 was masked
# 3 signif.x2   0.4 Effect of N+2 preview when N+1 was not masked
# 4 signif.x3   0.2 Effect of N+1 preview