Simulating observations for a 2-way ANOVA (in R) w/ mixed effects model and recovering true variance parameters [Gage R&R] I would like to create a simulation of a Gage R&R experiment in R. A Gage R&R is an experiment designed to analyze the variance contribution of several factors relative to the overall variance.  The context is often a measurement system where we'd like to know how much of the variation a measurement system is due to operator-to-operator variation, part-to-part variation, and random variation (repeatability) variation.  Observations from this type of experiments are typically modeled using a mixed effects model with a random effect for part, one for operator, a part:operator interaction, and a random error term. Note that each operator makes repeated measurements of the same part.

I'm trying to replicate the simulation described HERE where we specify the variance for each factor, generate observations, then fit a model and see how the estimates of the variance components compare with the true.  They show the general process but not the code or specifics for how to generate the data once the variances are specified.
if you already have the data, the process is pretty easy:
In R, the daewr package has a nice dataset to use as an example of fitting the model to existing data

library(lme4)
library(tidyverse)

#load data
data(gagerr)

#fit model
mod <- lmer(y ~ (1|part) + (1|oper) + (1|part:oper), data = gagerr)

#see variance of random effects
summary(mod)

Linear mixed model fit by REML ['lmerMod']
Formula: y ~ (1 | part) + (1 | oper) + (1 | part:oper)
   Data: gagerr

REML criterion at convergence: -133.9

Scaled residuals: 
     Min       1Q   Median       3Q      Max 
-2.43502 -0.36558 -0.01169  0.38978  1.94191 

Random effects:
 Groups    Name        Variance  Std.Dev.
 part:oper (Intercept) 0.0124651 0.11165 
 part      (Intercept) 0.0225515 0.15017 
 oper      (Intercept) 0.0000000 0.00000 
 Residual              0.0007517 0.02742 

Now I'd like to set the variance and simulate observations (then run the above analysis and compare to inputs).  My question is, how can I use the model to generate observations if all I care about are setting the variances?  In the reference article, they assume all the random effects are zero with variance sigma^2:  N(0, sigma^2).  i don't think it's as simple as just doing rnorm(60, 0, var^.5) and then adding the terms because of the interaction term.   The interaction term confuses me.  Do I need a bunch of matrix math to make sure the interaction aligns appropriately with the random effects such that when I run the analysis I can get a reasonable estimate of the true variance components?  Or is it more simple than that?
Thank you for any help you can provide.
 A: You are basically on the right track.

i don't think it's as simple as just doing rnorm(60, 0, var^.5) and then adding the terms because of the interaction term.

Correct, so you just have to simulate for the variance of the interaction too.
I find that the easiest way to simulate data for a mixed model do this is to use the model matrix, $Z$ for the random effect. Remember that the general equation for a mixed model is:
$$
Y = X\beta+Zb+e
$$
But here we have no fixed effects so it's just:
$$
Y = Zb+e
$$
where $Z$ is the model matrix the random effects and $b$ is the random effects coefficients vector
The problem is that unless the random structure is very simple, it can be quite tedious to construct $Z$ by hand. But, happily, there is an easy solution - just let the software do it for you. Here is an example using data corresponding to the model output in your question.
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

So here we just created the data frame for the factors, and added purely random noise to it to create a Y variable and used lFormula from the lme4 package to process the formula against the data without attempting to fit the model. During this processing the $Z$ model matrix is constructed and it's inverse $Zt$ is stored in the resulting object, so the last line there just transposes it to get $Z$.
Now we simulate the random effects themselves where I used standard deviations of 4, 3 and 2 for the random intercepts.
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)  

I had to check the order that these should go in. There are some rules for this in the documentation but I simply ran the code with 2 oper and 2 part and ran a full lmer model then extracted the random effects with ranef() and compared that the getME(mymodel, "b") which made it obvious. If this is confusing let me know and I'll add the code and output for that too.
Then we just simulate the outcome (with a unit level variance of 1) and fit the lmer model:
> 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   


And we see that we have nicely recovered the parameters 4, 3, 2, and 1 as the variance components
