library("MASS")
library("lme4")
set.seed(123)

simulate_data <- function(n_items = 60, n_subj = 27, n_subj_item_rep = 30,
                          mu = 308,
                          contrast1 = c(2/3, -1/3, -1/3), contrast2 = c(-1/3, 2/3, -1/3),
                          effect1 = 12, effect2 = -4.5,
                          subj_interc_sd = 43,
                          subj_slope1_sd = 17,
                          subj_slope2_sd = 112,
                          subj_cor_s1_i = 0.4,
                          subj_cor_s2_i = -0.65,
                          subj_cor_s1_s2 = -0.55,
                          item_interc_sd = 18,
                          item_slope1_sd = 6,
                          item_slope2_sd = 6,
                          item_cor_s1_i = 0.4,
                          item_cor_s2_i = 0.3,
                          item_cor_s1_s2 = -0.7,
                          subj_item_interc_sd = 20,
                          residual_sd = 180) {

  # generate design
  dat <- expand.grid(item = seq_len(n_items), subject = seq_len(n_subj))
  dat$item_group <- dat$item %% 3 + 1
  dat$subject_group <- dat$subject %% 3 + 1
  dat$condition <- (dat$item_group + dat$subject_group) %% 3 + 1

  # generate appropriate contrasts
  # R requires the generalized inverse of the contrast matrix
  contrast_mat <- ginv(rbind(contrast1, contrast2))
  dat$contrast1 <- contrast_mat[, 1][dat$condition]
  dat$contrast2 <- contrast_mat[, 2][dat$condition]

  # fixed grand mean
  dat$mean_response <- mu

  # condition effects
  dat$contrast_effect1 <- effect1 * dat$contrast1
  dat$contrast_effect2 <- effect2 * dat$contrast2

  # random effects subjects: intercept, slope1, slope2
  subject_effect <- mvrnorm(n_subj,
                            mu = c(0, 0, 0),
                            Sigma = matrix(c(subj_interc_sd^2,
                                             subj_cor_s1_i * subj_slope1_sd * subj_interc_sd,
                                             subj_cor_s2_i * subj_slope2_sd * subj_interc_sd,
                                             subj_cor_s1_i * subj_slope1_sd * subj_interc_sd,
                                             subj_slope1_sd^2,
                                             subj_cor_s1_s2 * subj_slope1_sd * subj_slope2_sd,
                                             subj_cor_s2_i * subj_slope2_sd * subj_interc_sd,
                                             subj_cor_s1_s2 * subj_slope1_sd * subj_slope2_sd,
                                             subj_slope2_sd^2), nrow = 3),
                            empirical = TRUE)  # mu & sigma as empirical (not population!) parameters for testing

  dat$subject_intercept <- subject_effect[dat$subject, 1]
  dat$subject_slope1 <- subject_effect[dat$subject, 2] * dat$contrast1
  dat$subject_slope2 <- subject_effect[dat$subject, 3] * dat$contrast2

  # random effects items: intercept, slope1, slope2
  item_effect <- mvrnorm(n_items,
                         mu = c(0, 0, 0),
                         Sigma = matrix(c(item_interc_sd^2,
                                          item_cor_s1_i * item_slope1_sd * item_interc_sd,
                                          item_cor_s2_i * item_slope2_sd * item_interc_sd,
                                          item_cor_s1_i * item_slope1_sd * item_interc_sd,
                                          item_slope1_sd^2,
                                          item_cor_s1_s2 * item_slope1_sd * item_slope2_sd,
                                          item_cor_s2_i * item_slope2_sd * item_interc_sd,
                                          item_cor_s1_s2 * item_slope1_sd * item_slope2_sd,
                                          item_slope2_sd^2), nrow = 3),
                         empirical = TRUE)

  dat$item_intercept <- item_effect[dat$item, 1]
  dat$item_slope1 <- item_effect[dat$item, 2] * dat$contrast1
  dat$item_slope2 <- item_effect[dat$item, 3] * dat$contrast2

  # random effects subject-item combinations
  dat$subj_item_intercept <- 0
  if (n_subj_item_rep > 1) {  # if more than 1 observations of the same subject-item combination
    dat$subj_item_intercept <- mvrnorm(n_subj * n_items, mu = 0, Sigma = subj_item_interc_sd^2, empirical = TRUE)[, 1]
  }

  # generate subject-item replications according to number of observations of the same subject-item combination
  dat <- do.call("rbind", replicate(n_subj_item_rep, dat, simplify = FALSE))

  # residual variation
  dat$residual_sd <- mvrnorm(nrow(dat), mu = 0, Sigma = residual_sd^2, empirical = TRUE)[, 1]

  # calculate response
  dat$response <- dat$mean_response + dat$contrast_effect1 + dat$contrast_effect2 +
    dat$item_intercept + dat$item_slope1 + dat$item_slope2 +
    dat$subject_intercept + dat$subject_slope1 + dat$subject_slope2 +
    dat$subj_item_intercept +
    dat$residual_sd

  # convert categorical variables to factors
  dat[, 1:5] <- lapply(dat[, 1:5], factor)

  # set specified contrasts
  contrasts(dat$condition) <- contrast_mat

  dat
}

d <- simulate_data()

# if more than 1 observations of the same subject-item combination ..
summary(lmer(response ~ contrast1 + contrast2 + (contrast1 + contrast2|item) +
               (contrast1 + contrast2|subject) + (1|subject:item), d))
# .. otherwise
# summary(lmer(response ~ contrast1 + contrast2 + (contrast1 + contrast2|item) +
#                (contrast1 + contrast2|subject), d))

Linear mixed model fit by REML ['lmerMod']
Formula: response ~ contrast1 + contrast2 + (contrast1 + contrast2 | item) +
    (contrast1 + contrast2 | subject) + (1 | subject:item)
   Data: d

REML criterion at convergence: 643678

Scaled residuals:
    Min      1Q  Median      3Q     Max
-4.1695 -0.6709 -0.0006  0.6728  4.2132

Random effects:
 Groups       Name        Variance Std.Dev. Corr
 subject:item (Intercept)   423.44  20.578
 item         (Intercept)   337.69  18.376
              contrast1      36.02   6.002   0.38
              contrast2      66.85   8.176   0.11 -0.88
 subject      (Intercept)  1819.74  42.658
              contrast1     270.89  16.459   0.45
              contrast2   12677.59 112.595  -0.65 -0.60
 Residual                 32431.73 180.088
Number of obs: 48600, groups:  subject:item, 1620; item, 60; subject, 27

Fixed effects:
            Estimate Std. Error t value
(Intercept)  308.000      8.600   35.82
contrast1     13.771      3.534    3.90
contrast2     -4.708     21.737   -0.22

Correlation of Fixed Effects:
          (Intr) cntrs1
contrast1  0.410
contrast2 -0.619 -0.557