Simular observações para um ANOVA (em R) de 2 vias com modelo de efeitos mistos e recuperar parâmetros de variância verdadeiros [Medidor R&R]

Aug 18 2020

Eu gostaria de criar uma simulação de um experimento Gage R&R em R. Um Gage R&R é um experimento projetado para analisar a contribuição da variância de vários fatores em relação à variância geral. O contexto é geralmente um sistema de medição onde gostaríamos de saber quanto da variação de um sistema de medição é devido à variação de operador para operador, variação de parte a parte e variação de variação aleatória (repetibilidade). As observações desse tipo de experimento são normalmente modeladas usando um modelo de efeitos mistos com um efeito aleatório para parte, um para o operador, uma parte: interação do operador e um termo de erro aleatório. Observe que cada operador faz medições repetidas da mesma peça.

Estou tentando replicar a simulação descrita AQUI, onde especificamos a variância para cada fator, geramos observações e, em seguida, ajustamos um modelo e vemos como as estimativas dos componentes da variância se comparam com o verdadeiro. Eles mostram o processo geral, mas não o código ou as especificações de como gerar os dados uma vez que as variações sejam especificadas.

se você já tem os dados, o processo é muito fácil:

Em R, o pacote daewr tem um bom conjunto de dados para usar como exemplo de ajuste do modelo aos dados 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

Agora, gostaria de definir a variância e simular observações (em seguida, execute a análise acima e compare as entradas). Minha pergunta é: como posso usar o modelo para gerar observações se tudo que me interessa é definir as variâncias? No artigo de referência, eles assumem que todos os efeitos aleatórios são zero com variância sigma ^ 2: N (0, sigma ^ 2). não acho que seja tão simples quanto apenas fazer rnorm (60, 0, var ^ .5) e adicionar os termos por causa do termo de interação. O termo de interação me confunde. Preciso de um monte de matemática de matriz para garantir que a interação se alinhe apropriadamente com os efeitos aleatórios de forma que, quando eu executar a análise, eu possa obter uma estimativa razoável dos verdadeiros componentes de variância? Ou é mais simples do que isso?

Obrigado por qualquer ajuda que você possa fornecer.

Respostas

5 RobertLong Aug 19 2020 at 13:02

Você está basicamente no caminho certo.

não acho que seja tão simples quanto apenas fazer rnorm (60, 0, var ^ .5) e adicionar os termos por causa do termo de interação.

Correto, então você só precisa simular a variação da interação também.

Acho que a maneira mais fácil de simular dados para um modelo misto é usar a matriz do modelo, $Z$para o efeito aleatório. Lembre-se de que a equação geral para um modelo misto é:

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

Mas aqui não temos efeitos fixos, então é apenas:

$$ Y = Zb+e $$

Onde $Z$ é a matriz do modelo, os efeitos aleatórios e $b$ é o vetor de coeficientes de efeitos aleatórios

O problema é que, a menos que a estrutura aleatória seja muito simples, pode ser entediante construir $Z$à mão. Mas, felizmente, existe uma solução fácil - deixe o software fazer isso por você. Aqui está um exemplo usando dados correspondentes à saída do modelo em sua pergunta.

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

Então, aqui acabamos de criar o quadro de dados para os fatores e adicionamos ruído puramente aleatório a ele para criar uma variável Y e usar lFormulado lme4pacote para processar a fórmula em relação aos dados sem tentar ajustar o modelo. Durante este processamento, a matriz do modelo $ Z $ é construída e seu inverso $ Zt $ é armazenado no objeto resultante, então a última linha apenas a transpõe para obter $ Z $ .

Agora simulamos os próprios efeitos aleatórios, onde usei desvios-padrão de 4, 3 e 2 para as interceptações aleatórias.

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)

Tive que verificar a ordem em que eles deveriam entrar. Existem algumas regras para isso na documentação, mas simplesmente executei o código com 2 opere 2 parte executei um lmermodelo completo , em seguida extraí os efeitos aleatórios ranef()e comparei o getME(mymodel, "b")que tornou óbvio . Se estiver confuso, me avise e adicionarei o código e a saída para isso também.