library(lme4)

MixLogFaker = function(N=1000,sat=F,nsub=5,nset=5,bias=0.5,nfact=2,ord=T,byit=F) {

# MixLogFaker simulates small-scale judgment experiments
# and estimates actual probability of getting p < .05

# Defaults
# N = 1000: Number of simulations
# sat = F: Test for satiation or not (interactions with Order)
# nsub = 5: Number of subjects (people making judgments)
# nset = 5: Number of sentence sets
# bias = 0.5: Higher means greater probability of 1 judgment
# nfact = 2: Number of factors (can only be 1 or 2)
# ord = T: Include Order as a covariate
# byit = F: Don't include Item as random variable

nitem = nset*nfact^2
totaln = nsub*nitem
ordran = numeric(nitem)
ordval = numeric(nitem)
if (ord && nfact == 1 && !sat) {nump = 2}
if (ord && nfact == 1 && sat) {nump = 3}
if (ord && nfact == 2 && !sat) {nump = 4}
if (ord && nfact == 2 && sat) {nump = 7}
if (!ord && nfact == 1) {nump = 1}
if (!ord && nfact == 2) {nump = 3}
sig = numeric(nump)
for (s in 1:N) {
 Judgment = as.numeric(runif(totaln)<bias)
 Subjects = NULL
 Items = NULL
 Factor1 = NULL
 Factor2 = NULL
 Sets = NULL
 Order = numeric(nitem)
 for (subj in 1:nsub) {
  if (ord) {
   # Create random ordering for this subject
   Order = NULL
   Order = c(Order,sample(1:nitem))
  }
  # Item numbers and factors
  if (nfact == 2) {
   item = 0
   for (set in 1:nset) {
    for (i in 0:1) {
     fact1 = i
     if (fact1 == 0) {fact1 = -1}
     for (j in 0:1) {
      fact2 = j
      if (fact2 == 0) {fact2 = -1}
      item = item + 1
      Items = c(Items,item)
      Subjects = c(Subjects,subj)
      Factor1 = c(Factor1,fact1)
      Factor2 = c(Factor2,fact2)
      Sets = c(Sets,set)
     }
    }
   }
  }
  if (nfact == 1) {
   item = 0
   for (set in 1:nset) {
    for (i in 0:1) {
     fact1 = i
     if (fact1 == 0) {fact1 = -1}
     item = item + 1
     Items = c(Items,item)
     Subjects = c(Subjects,subj)
     Factor1 = c(Factor1,fact1)
     Sets = c(Sets,set)
    }
   }
  }
 }
 if (!byit) {
  if (nfact == 1) {
   FakeData <- data.frame(Subjects,Items,Order,Factor1,Judgment)
   if (sat && ord) {
    GLMM.analysis = try(lmer(Judgment~Factor1*Order+(1|Subjects),
    data=FakeData,family="binomial"),T)
   }
   if (!sat && ord) {
    GLMM.analysis = try(lmer(Judgment~Factor1+Order+(1|Subjects),
    data=FakeData,family="binomial"),T)
   }
   if (!ord) {
    GLMM.analysis = try(lmer(Judgment~Factor1+(1|Subjects),
    data=FakeData,family="binomial"),T)
   }
  }
  if (nfact == 2) {
   FakeData = data.frame(Subjects,Items,Order,Factor1,Factor2,Judgment)
   if (sat && ord) {
    GLMM.analysis = try(lmer(Judgment~Factor1*Factor2*Order+(1|Subjects),
    data=FakeData,family="binomial"),T)
   }
   if (!sat && ord) {
    GLMM.analysis = try(lmer(Judgment~Factor1*Factor2+Order+(1|Subjects),
    data=FakeData,family="binomial"),T)
   }
   if (!ord) {
    GLMM.analysis = try(lmer(Judgment~Factor1*Factor2+(1|Subjects),
    data=FakeData,family="binomial"),T)
   }
  }
 }
 if (byit) {
  if (nfact == 1) {
   FakeData <- data.frame(Subjects,Items,Order,Factor1,Judgment)
   if (sat && ord) {
    GLMM.analysis = try(lmer(Judgment~Factor1*Order+(1|Subjects)+(1|Items),
    data=FakeData,family="binomial"),T)
   }
   if (!sat && ord) {
    GLMM.analysis = try(lmer(Judgment~Factor1+Order+(1|Subjects)+(1|Items),
    data=FakeData,family="binomial"),T)
   }
   if (!ord) {
    GLMM.analysis = try(lmer(Judgment~Factor1+(1|Subjects)+(1|Items),
    data=FakeData,family="binomial"),T)
   }
  }
  if (nfact == 2) {
   FakeData = data.frame(Subjects,Items,Order,Factor1,Factor2,Judgment)
   if (sat && ord) {
    GLMM.analysis = try(lmer(Judgment~Factor1*Factor2*Order+(1|Subjects)+(1|Items),
    data=FakeData,family="binomial"),T)
   }
   if (!sat && ord) {
    GLMM.analysis = try(lmer(Judgment~Factor1*Factor2+Order+(1|Subjects)+(1|Items),
    data=FakeData,family="binomial"),T)
   }
   if (!ord) {
    GLMM.analysis = try(lmer(Judgment~Factor1*Factor2+(1|Subjects)+(1|Items),
    data=FakeData,family="binomial"),T)
   }
  }
 }
 if (!is.character(try(summary(GLMM.analysis),T))) {
  values = summary(GLMM.analysis)@coefs
  for (f in 1:nump) {
   if (as.numeric(values[f+1,4]) < .05) {sig[f] = sig[f] + 1}
  }
 }
}
if (ord && nfact == 1 && !sat) {labs = c("Fact1","Order")}
if (ord && nfact == 1 && sat) {labs = c("Fact1","Order","Fact1xOrd")}
if (ord && nfact == 2 && !sat) {labs = c("Fact1","Fact2","Order","Fact1xFact2")}
if (ord && nfact == 2 && sat) {labs = c("Fact1","Fact2","Order","Fact1xFact2","Fact1xOrd","Fact2xOrd","F1xF2xOrd")}
if (!ord && nfact == 1) {labs = c("Fact1")}
if (!ord && nfact == 2) {labs = c("Fact1","Fact2","Fact1xFact2")}
alpha = sig/N
answer = rbind(labs,alpha)
return(answer)
}