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 MASS:mvrnom
function.
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
contrast.treatment
)
effect_A <- with(subset(d, condition == "A"), mean(b1 * X + noise))
all.equal(fixef(m1)[[1]], effect_A)
[1] TRUE
effect_B <- with(subset(d, condition == "B"), mean(b1 * X + noise))
all.equal(fixef(m1)[[2]], effect_B - effect_A)
[1] 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.