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?