Simulación de observaciones para un ANOVA de 2 vías (en R) con modelo de efectos mixtos y recuperación de parámetros de varianza real [Gage R&R]

Aug 18 2020

Me gustaría crear una simulación de un experimento R&R del sistema de medición en R. Un R&R del sistema de medición es un experimento diseñado para analizar la contribución de la varianza de varios factores en relación con la varianza general. El contexto es a menudo un sistema de medición en el que nos gustaría saber qué parte de la variación de un sistema de medición se debe a la variación de operador a operador, la variación de parte a parte y la variación de variación aleatoria (repetibilidad). Las observaciones de este tipo de experimentos se modelan típicamente utilizando un modelo de efectos mixtos con un efecto aleatorio para la parte, uno para el operador, una parte: interacción del operador y un término de error aleatorio. Tenga en cuenta que cada operador realiza mediciones repetidas de la misma pieza.

Estoy tratando de replicar la simulación descrita AQUÍ donde especificamos la varianza para cada factor, generamos observaciones, luego ajustamos un modelo y vemos cómo las estimaciones de los componentes de la varianza se comparan con la verdadera. Muestran el proceso general, pero no el código ni los detalles sobre cómo generar los datos una vez que se especifican las variaciones.

si ya tiene los datos, el proceso es bastante sencillo:

En R, el paquete daewr tiene un buen conjunto de datos para usar como ejemplo de cómo ajustar el modelo a los datos existentes.

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

Ahora me gustaría establecer la varianza y simular las observaciones (luego ejecutar el análisis anterior y compararlo con las entradas). Mi pregunta es, ¿cómo puedo usar el modelo para generar observaciones si lo único que me importa es establecer las variaciones? En el artículo de referencia, asumen que todos los efectos aleatorios son cero con varianza sigma ^ 2: N (0, sigma ^ 2). No creo que sea tan simple como hacer rnorm (60, 0, var ^ .5) y luego agregar los términos debido al término de interacción. El término de interacción me confunde. ¿Necesito un montón de matemáticas matriciales para asegurarme de que la interacción se alinea adecuadamente con los efectos aleatorios de modo que cuando ejecuto el análisis pueda obtener una estimación razonable de los componentes de la varianza real? ¿O es más simple que eso?

Gracias por cualquier ayuda que pueda brindar.

Respuestas

5 RobertLong Aug 19 2020 at 13:02

Básicamente estás en el camino correcto.

No creo que sea tan simple como hacer rnorm (60, 0, var ^ .5) y luego agregar los términos debido al término de interacción.

Correcto, así que solo tienes que simular la varianza de la interacción también.

Encuentro que la forma más fácil de simular datos para un modelo mixto es usar la matriz del modelo, $Z$para el efecto aleatorio. Recuerde que la ecuación general para un modelo mixto es:

$$ Y = X\beta+Zb+e $$

Pero aquí no tenemos efectos fijos, por lo que es solo:

$$ Y = Zb+e $$

dónde $Z$ es la matriz del modelo los efectos aleatorios y $b$ es el vector de coeficientes de efectos aleatorios

El problema es que, a menos que la estructura aleatoria sea muy simple, puede resultar bastante tedioso de construir $Z$manualmente. Pero, afortunadamente, hay una solución fácil: deje que el software lo haga por usted. Aquí hay un ejemplo que usa datos correspondientes a la salida del modelo en su pregunta.

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

Entonces, aquí solo creamos el marco de datos para los factores y le agregamos ruido puramente aleatorio para crear una variable Y y lo usamos lFormuladel lme4paquete para procesar la fórmula contra los datos sin intentar ajustar el modelo. Durante este procesamiento , se construye la matriz del modelo $ Z $ y su inverso $ Zt $ se almacena en el objeto resultante, por lo que la última línea allí simplemente lo transpone para obtener $ Z $ .

Ahora simulamos los efectos aleatorios en sí mismos donde usé desviaciones estándar de 4, 3 y 2 para las intersecciones aleatorias.

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)

Tuve que verificar el orden en que deberían entrar. Hay algunas reglas para esto en la documentación, pero simplemente ejecuté el código con 2 opery 2 party ejecuté un lmermodelo completo , luego extraje los efectos aleatorios ranef()y comparé los getME(mymodel, "b")que lo hicieron obvio. . Si esto es confuso, avíseme y agregaré el código y la salida para eso también.