Simulating Outcomes from a Random Effects Model - Strange Behavior when Changing Number of Effects

Question

I've have been doing my best to generate observations from a random effects model so that I can compare estimates of parameters to true parameters for a variety of conditions (like number of random effects, different magnitudes of standard deviation among the effects, etc).

I have been following the guidance from this post (my last post) where I got some good advice on how to go about simulating from these types of models. In summary, the advice was to set up a dummy experiment, extract the design matrix Z from that setup, draw your random effects and assemble into a vector, then use the design matrix Z and the random effect to construct the simulated observations. Then fit a model to the observations to see the estimate:

set.seed(15)
n.part <- 20  # number of parts
n.oper <- 20  # number of opers
n.reps <- 2   # number of replications

dt <- expand.grid(part = LETTERS[1:n.part], oper = 1:n.oper, reps = 1:n.reps)

dt$Y <- 10 + rnorm(n.part*n.oper*n.reps)

myformula <- "Y ~ (1|part) + (1|oper) + (1|part:oper)"  # model formula

mylF <- lFormula(eval(myformula), data = dt) # Process the formula against the data
Z <- mylF$reTrms$Zt %>% as.matrix() %>% t()  # Extract the Z matrix

b1 <- rnorm(n.part * n.oper, 0 , 4)   # random interecepts for the interaction
b2 <- rnorm(n.oper, 0, 3)             # random interecepts for oper
b3 <- rnorm(n.part, 0, 2)             # random interecepts for part

b <- c(b1, b2, b3)  

dt$Y <- 10 + Z %*% b + rnorm(nrow(dt))
> lmer(eval(myformula), data = dt ) %>% summary()
Linear mixed model fit by REML ['lmerMod']
Formula: Y ~ (1 | part) + (1 | oper) + (1 | part:oper)
   Data: dt

REML criterion at convergence: 3776.8

Scaled residuals: 
     Min       1Q   Median       3Q      Max 
-2.42747 -0.46098  0.01696  0.46941  2.44928 

Random effects:
 Groups    Name        Variance Std.Dev.
 part:oper (Intercept) 16.833   4.103   
 oper      (Intercept) 10.183   3.191   
 part      (Intercept)  4.840   2.200   
 Residual               1.009   1.005   

I have now been running simulations where I hold the st_dev of the random effects for: oper, and part:oper constant and vary the magnitude of the effect of part. I am seeing some behavior that I do not understand: If I use a equal number of parts and operators, for example 10 and 10, then I can recover the true parameters for standard deviation across a wide range of sd for part. However, if I change the number of part and operators to, for example 10 and 9, then the results get very wonky and I cannot recover the right params for sd of part or operator. One misses and one low. This does not appear to be an affect of just "sample size"...if i increase the number of both parts and operators but make them slightly different from one another I still see this same behavior (example: 20 parts, 19 oper)

See the following images: this first is a simulated experiment with n=10 parts, n=10 oper. The red dots are the true population standard deviations for those effects.

This 2nd is n=10 parts, n=9 oper. Again, the red dots are true pop parameters.

20 parts, 19 oper:

Is this to be expected for designs as I described? Or is there likely an error in the code for my simulations? Perhaps I can't just extract a design matrix so simply as described in the previous post?

Answer 1

The problem here seems to be that the order of the random effects in the Z matrix will not always be the same, so you can't always rely on the line b <- c(b1, b2, b3) being correct. It appears that lme4 constructs the model matrix so that it expects the vector of ranom effects b to be in descending order, from the one with the most levels to the one with the fewest. In your case, you have the part:oper interactions as a grouping variable so this will always have the most levels, so b1 should come first. When there are equal numbers, such as in the code block above it appears that lme4 might be using alphabetical order, but when you reduce the number forpart to below that of oper then part will come next and you will need to use b <- c(b1, b3, b2). So something like this should do the trick:

if (n.part < n.oper) {
  b <- c(b1, b3, b2)  
} else {
  b <- c(b1, b2, b3)  
}