Simulazione di osservazioni per un ANOVA a 2 vie (in R) con modello a effetti misti e recupero dei parametri di varianza reali [Gage R&R]

Aug 18 2020

Vorrei creare una simulazione di un esperimento di Gage R&R in R. Un Gage R&R è un esperimento progettato per analizzare il contributo della varianza di diversi fattori rispetto alla varianza complessiva. Il contesto è spesso un sistema di misurazione in cui vorremmo sapere quanta variazione di un sistema di misurazione è dovuta alla variazione da operatore a operatore, da parte a parte e da variazione casuale (ripetibilità). Le osservazioni di questo tipo di esperimenti sono tipicamente modellate utilizzando un modello a effetti misti con un effetto casuale per parte, uno per operatore, una parte: interazione operatore e un termine di errore casuale. Notare che ogni operatore effettua misurazioni ripetute della stessa parte.

Sto cercando di replicare la simulazione descritta QUI in cui specifichiamo la varianza per ciascun fattore, generiamo osservazioni, quindi adattiamo un modello e vediamo come le stime delle componenti della varianza si confrontano con il vero. Mostrano il processo generale ma non il codice o le specifiche su come generare i dati una volta specificate le varianze.

se hai già i dati, il processo è abbastanza semplice:

In R, il pacchetto daewr ha un bel dataset da usare come esempio per adattare il modello ai dati esistenti

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

Ora vorrei impostare la varianza e simulare le osservazioni (quindi eseguire l'analisi sopra e confrontare con gli input). La mia domanda è: come posso utilizzare il modello per generare osservazioni se tutto ciò che mi interessa è impostare le varianze? Nell'articolo di riferimento, assumono che tutti gli effetti casuali siano zero con varianza sigma ^ 2: N (0, sigma ^ 2). Non penso sia semplice come fare semplicemente rnorm (60, 0, var ^ .5) e poi aggiungere i termini a causa del termine di interazione. Il termine di interazione mi confonde. Ho bisogno di un po 'di matematica a matrice per assicurarmi che l'interazione si allinei in modo appropriato con gli effetti casuali in modo tale che quando eseguo l'analisi posso ottenere una stima ragionevole delle componenti della varianza reale? O è più semplice di così?

Grazie per tutto l'aiuto che puoi fornire.

Risposte

5 RobertLong Aug 19 2020 at 13:02

Sei fondamentalmente sulla strada giusta.

Non penso sia semplice come fare semplicemente rnorm (60, 0, var ^ .5) e poi aggiungere i termini a causa del termine di interazione.

Esatto, quindi devi solo simulare anche la varianza dell'interazione.

Trovo che il modo più semplice per simulare i dati per un modello misto è usare la matrice del modello, $Z$per l'effetto casuale. Ricorda che l'equazione generale per un modello misto è:

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

Ma qui non abbiamo effetti fissi quindi è solo:

$$ Y = Zb+e $$

dove $Z$ è la matrice del modello gli effetti casuali e $b$ è il vettore dei coefficienti degli effetti casuali

Il problema è che, a meno che la struttura casuale non sia molto semplice, può essere piuttosto noiosa da costruire $Z$a mano. Ma, fortunatamente, esiste una soluzione semplice: lascia che sia il software a farlo per te. Ecco un esempio che utilizza i dati corrispondenti all'output del modello nella tua domanda.

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

Quindi qui abbiamo appena creato il data frame per i fattori e aggiunto un rumore puramente casuale per creare una variabile Y e utilizzato lFormuladal lme4pacchetto per elaborare la formula rispetto ai dati senza tentare di adattare il modello. Durante questa elaborazione viene costruita la matrice del modello $ Z $ e la sua inversa $ Zt $ viene memorizzata nell'oggetto risultante, quindi l'ultima riga la traspone per ottenere $ Z $ .

Ora simuliamo gli effetti casuali stessi dove ho usato deviazioni standard di 4, 3 e 2 per le intercettazioni casuali.

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)

Ho dovuto controllare l'ordine in cui questi dovrebbero entrare. Ci sono alcune regole per questo nella documentazione ma ho semplicemente eseguito il codice con 2 opere 2 parte ho eseguito un lmermodello completo , quindi ho estratto gli effetti casuali ranef()e getME(mymodel, "b")ho confrontato quello che lo rendeva ovvio . Se questo crea confusione fammelo sapere e aggiungerò il codice e l'output anche per quello.