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
