require(lme4)
require(plyr)

participant <- rep(1:12, 3)
essay <- c(rep(1, 9), rep(2, 9), rep(3, 9), rep(4, 9))
quality <- ifelse(essay == 1 | essay == 2, "high", "low")
race <- c("white", "black", "black", "white", "white", "white", "white", "white", "white",
          "white", "white", "white", "white", "white", "white", "white", "black", "black",
          "white", "black", "black", "white", "white", "white", "white", "white", "white",
          "white", "white", "white", "white", "white", "white", "white", "black", "black")
gender <- c("female", "male", "female", "male", "male", "male", "male", "male", "male",
            "male", "male", "male", "male", "male", "male", "female", "male", "female",
            "female", "male", "female", "male", "male", "male", "male", "male", "male",
            "male", "male", "male", "male", "male", "male", "female", "male", "female")

d <- data.frame(participant, essay, quality, race, gender)

for(i in 1:35)
{
  participant <- participant + 12
  essay <- essay + 4
  newdat <- data.frame(participant, essay, quality, race, gender)

  d <- rbind(d, newdat)
}

check_var <- function(var)
{
  tab <- table(var)
  newvar <- character()

  for(i in var)
  {
    if(i == names(tab[tab == 1]))
    {
      newvar <- c(newvar, "different")
    } else
    {
      newvar <- c(newvar, "same")
    }
  }

  return(newvar)
}

# Mark, for a given participant, which essay is "different"
d <- ddply(d, "participant", mutate, different = check_var(quality))

# Make each variable numeric for the purposes of the simulation
d$quality <- ifelse(d$quality == "low", -.5, .5)
d$race <- ifelse(d$race == "black", -.5, .5)
d$gender <- ifelse(d$gender == "female", -.5, .5)
d$different <- ifelse(d$different == "same", -.5, .5)

# Random seed
set.seed(2352)

# Number of simulations
reps <- 50
# Create a storage space for the effects
effs <- matrix(NA, ncol = 10, nrow = reps)

# For each simulation
for(i in 1:reps)
{
  # Fixed effects.  A quality effect and a contrast effect for quality
  d$score <- .5 * d$quality + 1 * d$different * d$quality

  # Random effects at the participant level
  d <- ddply(d, "participant", mutate, r_int = rnorm(1, sd = .5),
             r_q = rnorm(1, sd = .5),
             score = score + r_int + r_q * quality)

  # Random effects at the essay level
  d <- ddply(d, "essay", mutate, g_int = rnorm(1, sd = .5),
             g_r = rnorm(1, sd = .5),
             g_g = rnorm(1, sd = .5),
             g_r_g = rnorm(1, sd = .5),
             score = score + g_int + g_r * race + g_g * gender + g_r_g * race * gender)

  # Observation-level error
  d$score <- d$score + rnorm(dim(d)[1], sd = 1)

  # Fit the model
  mod <- lmer(score ~ race * gender * quality + different * quality + (race * gender | essay) + (quality | participant), data = d)

  # Store the coefficients
  colnames(effs) <- names(fixef(mod))
  effs[i, ] <- fixef(mod)

  # Print the current simulation
  print(i)
}

# Results
round(colMeans(effs), digits = 2)

        (Intercept)                race              gender             quality
               0.00               -0.03                0.02                0.50
          different         race:gender        race:quality      gender:quality
               0.01               -0.03                0.00                0.03
  quality:different race:gender:quality
               0.97               -0.02