# 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)

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?

• It's not clear what you want to do. Do you know the fixed effects and just want to simulate the outcome with different variances for the random effects ? Aug 18, 2020 at 20:33
• I would like to assign values to whatever parameters are necessary in order to generate simulated observations. Then I would like to treat these simulated observations as a random sample and fit a model and compare to the parameter estimates I just set. I will look at all the parameters but the ones I really care about are the variances of the random effects. The part that is confusing me is how to actually create the observations. Thank you. Aug 18, 2020 at 21:37
• Have you thought to simulate data using Stan? Like, make_stancode( y ~ (1|part) + (1|oper) + (1|part:oper), data = gagerr) and then use the produced Stan code for simulation. Dec 15, 2020 at 9:30

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

• Amazing answer - thank you so much. This is complicated as I expected but not so bad that i cant wrap my head around it with some study. Really appreciate such a thorough answer. Aug 21, 2020 at 2:44
• No problem, happy to help. The key is to understanding the general idea with design matrices. A good place to start is by looking at a simple linear model $Y = X\beta$ and then look a model with just random intercepts. With random effects the design matrix $Z$is more complicated but you just need to realise that all it does is map the random effect(s) for each component (each level of the grouping factor(s)) to the response. Also it will help if you understand how to multiply a matrix and a vector (eg a 2x2 matrix and a 2x1 vector) Aug 21, 2020 at 7:36
• sorry for the ping, but do you know if there is a way to use this same method when the number of effects is different in each random variable? I would think extracting the Z should still work, but I find if I change n.oper in the above example from 20 to 19 then this method doesn't seem to work - problem is exacerbated as you increase the magnitude of the sd for the random intercepts for part. Sep 3, 2020 at 4:02
• What do you mean "doesn't seem to work" ? Sep 3, 2020 at 5:49