# Simulating linear mixed model data for factorial design with 3 levels

I struggle to simulate linear mixed model data for the following within-subject within-item factorial design (with crossed participant and item effects) where the categorical independent variable has three levels:

\begin{aligned} &Y_{c,s,i}\ =\ \beta_0\ +\ S_{0,s}\ +\ I_{0,i}\ +\ C_{0,s,i}\ +\ \left(\beta_1\ +\ S_{1,s}\ +\ I_{1,i}\right)X_{1,c}\ +\left(\beta_2\ +\ S_{2,s}\ +\ I_{2,i}\right)X_{2,c}\ +\ \epsilon_{c,s,i},\\ &\begin{pmatrix} S_{0,s}\\ S_{1,s}\\ S_{2,s} \end{pmatrix} \sim N \begin{bmatrix} \begin{pmatrix} 0\\ 0\\ 0 \end{pmatrix}, \begin{pmatrix} \tau_{00}^2 & \rho_{s_{01}}\tau_{00}\tau_{11} & \rho_{s_{02}}\tau_{00}\tau_{22}\\ \rho_{s_{01}}\tau_{00}\tau_{11} & \tau_{11}^2 & \rho_{s_{12}}\tau_{11}\tau_{22}\\ \rho_{s_{02}}\tau_{00}\tau_{22} & \rho_{s_{12}}\tau_{11}\tau_{22} & \tau_{22}^2 \end{pmatrix} \end{bmatrix},\\ &\begin{pmatrix} I_{0,i}\\ I_{1,i}\\ I_{2,i} \end{pmatrix} \sim N \begin{bmatrix} \begin{pmatrix} 0\\ 0\\ 0 \end{pmatrix}, \begin{pmatrix} \omega_{00}^2 & \rho_{i_{01}}\omega_{00}\omega_{11} & \rho_{i_{02}}\omega_{00}\omega_{22}\\ \rho_{i_{01}}\omega_{00}\omega_{11} & \omega_{11}^2 & \rho_{i_{12}}\omega_{11}\omega_{22}\\ \rho_{i_{02}}\omega_{00}\omega_{22} & \rho_{i_{12}}\omega_{11}\omega_{22} & \omega_{22}^2 \end{pmatrix} \end{bmatrix},\\ &C_{0,s,i}\ \sim\ N \left(0,\ \phi^2\right),\\ &\epsilon_{c,s,i}\ \sim\ N \left(0,\ \sigma^2\right) \end{aligned}

$Y_{c,s,i}$ refers to the dependent variable, the subscripts stand for condition (c), subject (s) and item (i), the by-subject, by-item and by-subject-item-combination intercepts are $S_{0,s}$, $I_{0,i}$ and $C_{0,s,i}$. The corresponding by-subject and by-item slopes are $S_{1,s}$, $S_{2,s}$ and $I_{1,i}$, $I_{2,i}$, respectively. The fixed effects are $\beta_0$, $\beta_1$ and $\beta_2$ and the residual error is $\epsilon_{c,s,i}$.

I hope the mapping between this model formulation and simulata_data() is clear:

• $\beta_1$ -> effect1
• $\tau_{00}$ -> subj_interc_sd
• $\rho_{s_{12}}$ -> subj_cor_s1_s2
• $\omega_{11}$ -> item_slope1_sd
• $\phi$ -> subj_item_interc_sd
• ...

The default parameters I used correspond to the estimates I got from an analysis of previous data.
As I am using the mvrnorm() function with empirical = TRUE I expected that the estimated parameters from lmer and my specified parameters for the simulation would match closely.
However, they differ more than I expected:

• effect1: should be 12, estimated as 13.771 (Fixed effects table: contrast1 in summary below)
• item_cor_s2_i: should be 0.3, estimated as 0.11 (Random effects table: Corr(item, contrast2) in summary below))
• item_cor_s1_s2: should be -0.7, estimated as -0.88 (Random effects table: Corr(contrast1, contrast2) in summary below)

Since this is my first time programming a simulation for mixed models, I'm not sure if the (for me) unexpected results are caused by a programming error or if the default parameters are difficult to estimate for lmer.
I would be grateful if somebody more experienced in mixed model simulations could comment on this.

Running this code ...

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))


... I get the following output

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

• What is your question? – AdamO Mar 20 '18 at 15:18
• Thank you for providing runnable code. Unfortunately, I don't see the output you refer to in the R output. Which estimates are you unhappy about? That said, have you tried running the simulation with different seeds and looking how much the estimates vary? – Stephan Kolassa Mar 21 '18 at 11:23
• @Stephan Kolassa: I hope my edit makes this clear. The estimates do indeed vary. However, I expected that they would match closer with my specified parameters for the simulation. – statmerkur Mar 22 '18 at 16:57

Here is my attempt to simulate the data. I've tried to simplify it.

library(MASS)
library(mvtnorm)
library(dplyr)
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) {
contrast_mat <- ginv(rbind(contrast1, contrast2))
m_s <- diag(c(subj_interc_sd, subj_slope1_sd, subj_slope2_sd))
m_s[lower.tri(m_s)] <- c(subj_cor_s1_i, subj_cor_s2_i, subj_cor_s1_s2)
m_s[upper.tri(m_s)] <- c(subj_cor_s1_i, subj_cor_s2_i, subj_cor_s1_s2)
subject_effects <- rmvnorm(n_subj, sigma = sdcor2cov(m_s))
m_i <- diag(c(item_interc_sd, item_slope1_sd, item_slope2_sd))
m_i[lower.tri(m_i)] <- c(item_cor_s1_i, item_cor_s2_i, item_cor_s1_s2)
m_i[upper.tri(m_i)] <- c(item_cor_s1_i, item_cor_s2_i, item_cor_s1_s2)
item_effects <- rmvnorm(n_items, sigma = sdcor2cov(m_i))

expand.grid(
item = seq_len(n_items),
subject = seq_len(n_subj),
replicate = seq_len(n_subj_item_rep)
) %>%
mutate(
item_group = item %% 3 + 1,
subject_group = subject %% 3 + 1,
condition = (item_group + subject_group) %% 3 + 1,
contrast1 = contrast_mat[, 1][condition],
contrast2 = contrast_mat[, 2][condition],
fixed = mu + effect1 * contrast1 + effect2 * contrast2,
random_subject = subject_effects[subject, 1] +
subject_effects[subject, 2] * contrast1 +
subject_effects[subject, 3] * contrast2,
random_item = item_effects[item, 1] +
item_effects[item, 2] * contrast1 +
item_effects[item, 3] * contrast2,
item_subject = interaction(item, subject),
random_interaction = rnorm(n_items * n_subj, sd = subj_item_interc_sd)[item_subject],
noise = rnorm(n(), sd = residual_sd),
response = fixed + random_subject + random_item + random_interaction + noise
) %>%
mutate_at(c("item_group", "subject_group", "item", "subject", "condition"), factor) -> dat
contrasts(dat$condition) <- contrast_mat return(dat) }  Next we can simulate the data several times and model these datasets. In the example below I only extract the fixed effect parameters. set.seed(123) replicate( 10, { simulate_data() %>% lmer( formula = response ~ contrast1 + contrast2 + (contrast1 + contrast2|item) + (contrast1 + contrast2|subject) + (1|subject:item) ) %>% fixef() }) %>% t() -> sim summary(sim)  The output of the summary (below) indicates quite strong differences in estimates among the simulations. This isn't surprising given that several standard deviations are very high ($\sigma$= 180,$\tau_{22} = 112\$). As a result the signal-to-noise ratio is quite low.

  (Intercept)      contrast1       contrast2
Min.   :306.9   Min.   :12.32   Min.   :-25.122
1st Qu.:309.4   1st Qu.:13.59   1st Qu.:-19.757
Median :310.3   Median :16.05   Median : -7.593
Mean   :311.3   Mean   :15.70   Mean   : -6.242
3rd Qu.:310.6   3rd Qu.:17.61   3rd Qu.:  2.968
Max.   :320.1   Max.   :18.36   Max.   : 18.226

• Thanks for simplifying the code. AFAICS, there is no difference in how you calculate the response variable compared to my approach. So, do you thinks my coding is, apart from being overly complicated, correct? I explicitly used MASS::mvrnorm() with empirical = TRUE to specify empirical not population parameters for the multivariate normal distributions. As a result my estimates vary not as strongly as yours but still quite strongly. I think this is a result of the low signal-to-noise ratio you mentioned. Can you add these two issues to your answer? – statmerkur Mar 31 '18 at 11:59