Simulation d'observations pour un modèle ANOVA bidirectionnel (en R) avec effets mixtes et récupération des paramètres de variance réels [Gage R&R]

Aug 18 2020

Je voudrais créer une simulation d'une expérience Gage R&R en R. Un Gage R&R est une expérience conçue pour analyser la contribution à la variance de plusieurs facteurs par rapport à la variance globale. Le contexte est souvent un système de mesure où nous aimerions savoir quelle part de la variation un système de mesure est due à la variation d'opérateur à opérateur, à la variation de pièce à pièce et à la variation de variation aléatoire (répétabilité). Les observations de ce type d'expériences sont généralement modélisées à l'aide d'un modèle à effets mixtes avec un effet aléatoire pour la partie, un pour l'opérateur, une partie: interaction de l'opérateur et un terme d'erreur aléatoire. Notez que chaque opérateur effectue des mesures répétées de la même pièce.

J'essaie de reproduire la simulation décrite ICI où nous spécifions la variance pour chaque facteur, générons des observations, puis ajustons un modèle et voyons comment les estimations des composantes de la variance se comparent au vrai. Ils montrent le processus général mais pas le code ou les détails sur la façon de générer les données une fois les écarts spécifiés.

si vous avez déjà les données, le processus est assez simple:

Dans R, le package daewr a un joli jeu de données à utiliser comme exemple d'ajustement du modèle aux données existantes

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

Maintenant, j'aimerais définir la variance et simuler les observations (puis exécuter l'analyse ci-dessus et comparer aux entrées). Ma question est la suivante: comment puis-je utiliser le modèle pour générer des observations si tout ce qui m'intéresse est de définir les écarts? Dans l'article de référence, ils supposent que tous les effets aléatoires sont nuls avec une variance sigma ^ 2: N (0, sigma ^ 2). Je ne pense pas que ce soit aussi simple que de faire rnorm (60, 0, var ^ .5) et d'ajouter ensuite les termes à cause du terme d'interaction. Le terme d'interaction me déroute. Ai-je besoin d'un tas de mathématiques matricielles pour m'assurer que l'interaction s'aligne correctement avec les effets aléatoires de sorte que lorsque j'exécute l'analyse, je puisse obtenir une estimation raisonnable des composantes de la variance réelle? Ou est-ce plus simple que ça?

Merci pour toute aide que vous pouvez fournir.

Réponses

5 RobertLong Aug 19 2020 at 13:02

Vous êtes fondamentalement sur la bonne voie.

Je ne pense pas que ce soit aussi simple que de faire rnorm (60, 0, var ^ .5) et d'ajouter ensuite les termes à cause du terme d'interaction.

Correct, il vous suffit donc de simuler la variance de l'interaction.

Je trouve que le moyen le plus simple de simuler des données pour un modèle mixte consiste à utiliser la matrice du modèle, $Z$pour l'effet aléatoire. N'oubliez pas que l'équation générale pour un modèle mixte est:

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

Mais ici, nous n'avons pas d'effets fixes donc c'est juste:

$$ Y = Zb+e $$

où $Z$ est la matrice du modèle les effets aléatoires et $b$ est le vecteur des coefficients d'effets aléatoires

Le problème est qu'à moins que la structure aléatoire ne soit très simple, il peut être assez fastidieux de construire $Z$par la main. Mais, heureusement, il existe une solution simple: laissez simplement le logiciel le faire pour vous. Voici un exemple utilisant des données correspondant à la sortie du modèle dans votre 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

Donc, ici, nous avons juste créé la trame de données pour les facteurs, et y avons ajouté un bruit purement aléatoire pour créer une variable Y et utilisé à lFormulapartir du lme4package pour traiter la formule par rapport aux données sans essayer d'ajuster le modèle. Pendant ce traitement, la matrice du modèle $ Z $ est construite et son inverse $ Zt $ est stocké dans l'objet résultant, donc la dernière ligne la transpose juste pour obtenir $ Z $ .

Maintenant, nous simulons les effets aléatoires eux-mêmes où j'ai utilisé des écarts types de 4, 3 et 2 pour les interceptions aléatoires.

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)

J'ai dû vérifier l'ordre dans lequel ils devaient entrer. Il y a quelques règles pour cela dans la documentation, mais j'ai simplement exécuté le code avec 2 operet 2 partet j'ai exécuté un lmermodèle complet , puis j'ai extrait les effets aléatoires ranef()et comparé ceux qui le getME(mymodel, "b")rendaient évident . Si cela est déroutant, faites-le moi savoir et j'ajouterai le code et la sortie pour cela aussi.