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.