How to simulate "empirical" fixed effects in linear mixed models

Question

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.

Answer 1

Since the mean of the simulated residuals is not exactly zero, the model cannot possibly disentangle the mean residual from the true effect. Nevertheless, the model clearly does a good job in taking sampling variability into account, because both your estimates are less than one standard error from their true values. And if you increase the number of observations you can get as close as you want to the original values used in the simulation.

Edit: In your code the mean of the residual is zero but only when summed over the whole dataset, i.e. pooling together the two conditions A and B. The residuals' means are not zero when averaged separately in the two conditions, hence the difference between the estimates and the true generative values. I don't think it is possible to obtain the same exact values of the fixed-effects used in the simulation, unless you force the residuals to sum to zero within conditions, e.g. code below.

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]

# simulate residuals separately to make sure they sum to zero within conditions
d$noise[d$condition=="A"] <- MASS::mvrnorm(nrow(d[d$condition=="A",]), 
                         mu = 0, 
                         Sigma = 3^2*(60-1)/(60-2), 
                         empirical = TRUE)

d$noise[d$condition=="B"] <- MASS::mvrnorm(nrow(d[d$condition=="B",]), 
                         mu = 0, 
                         Sigma = 3^2*(60-1)/(60-2), 
                         empirical = TRUE)

d <- within(d, y <- b1 * X + subject_intercept + noise)

# now the fixed-effect estimates are the same as the true values, up to some numerical error
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: 327.9733
Random effects:
 Groups   Name        Std.Dev.
 subject  (Intercept) 4.635   
 Residual             3.141   
Number of obs: 60, groups:  subject, 10
Fixed Effects:
(Intercept)   conditionB  
  -3.44e-14     5.00e+01 

Edit 2: When setting empirical=TRUE the values generated are such that the sample variance computed using Bessel correction (that is by using of $n-1$ instead of $n$ in the denominator) correspond exactly to the one used for generating the values. In the code above the residuals are generated for two smaller vectors of length 30 instead of a single one of length 60, hence the estimated sample variance would be different from the intended value, unless this value is corrected as pointed out by @amoeba in his comment. In this particular case one needs to multiply the desired value by $\frac{60-1}{60-2}$, and more generally if one wants to generate a sample of length $n$ by concatenating $k$ vectors of equal length, the standard deviation used to generate each of the $k$ vectors should be multiplied by the factor $\frac{n-1}{n-k}$.