Is there a way to simulate a linear mixed model where the estimated fixed effects match the fixed effects one specifies in the simulation?
At first I assumed that the following code would achieve this.
Please note: I set
empirical = TRUE within the
library("lme4") set.seed(123) d <- data.frame("subject" = gl(10, 6), "condition" = gl(2, 1, 60, labels = c("A", "B")), "X" = rep(0:1, 30), "b1" = 50) d$subject_intercept <- MASS::mvrnorm(nlevels(d$subject), mu = 0, Sigma = 5^2, empirical = TRUE) [d$subject] d$noise <- MASS::mvrnorm(nrow(d), mu = 0, Sigma = 3^2, empirical = TRUE) d <- within(d, y <- b1 * X + subject_intercept + noise) print(m1 <- lmer(y ~ condition + (1|subject), d)) Linear mixed model fit by REML ['lmerMod'] Formula: y ~ condition + (1 | subject) Data: d REML criterion at convergence: 323.7271 Random effects: Groups Name Std.Dev. subject (Intercept) 4.412 Residual 3.035 Number of obs: 60, groups: subject, 10 Fixed Effects: (Intercept) conditionB -0.2248 50.4496
However, the fixed effect estimates for the intercept and the contrast ("conditionB") are:
- intercept: the simulated (constant) effect of condition A + the mean of the simulated residuals within this condition
- contrast: the simulated (constant) effect of condition B + the mean of the simulated residuals within this condition - the intercept (because of
effect_A <- with(subset(d, condition == "A"), mean(b1 * X + noise)) all.equal(fixef(m1)[], effect_A)  TRUE effect_B <- with(subset(d, condition == "B"), mean(b1 * X + noise)) all.equal(fixef(m1)[], effect_B - effect_A)  TRUE
This seems to be a problem because AFAIK one usually only specifies the residual variation and the mean of the residuals (0) across conditions.