Karma efektli 2 yönlü ANOVA (R'de) için gözlemleri simüle etme ve gerçek varyans parametrelerini kurtarma [Gage R&R]

Question

R'de bir Gage Ar-Ge deneyinin bir simülasyonunu oluşturmak istiyorum. Bir Gage Ar-Ge, genel varyansa göre çeşitli faktörlerin varyans katkısını analiz etmek için tasarlanmış bir deneydir. Bağlam, genellikle bir ölçüm sisteminin varyasyonunun ne kadarının operatörden operatöre, parçadan kısma varyasyona ve rastgele varyasyon (tekrarlanabilirlik) varyasyonuna bağlı olduğunu bilmek istediğimiz bir ölçüm sistemidir. Bu tür deneylerden elde edilen gözlemler, tipik olarak parça için rastgele efekt, biri operatör için, bir parça: operatör etkileşimi ve rastgele hata terimiyle karışık etkiler modeli kullanılarak modellenir. Her operatörün aynı parça üzerinde tekrarlayan ölçümler yaptığını unutmayın.

BURADA açıklanan simülasyonu, her faktör için varyansı belirlediğimiz, gözlemler ürettiğimiz, ardından bir model uydurduğumuz ve varyans bileşenlerinin tahminlerinin doğru ile nasıl karşılaştırıldığını gördüğümüz simülasyonu kopyalamaya çalışıyorum . Genel süreci gösterirler, ancak varyanslar belirlendikten sonra verilerin nasıl üretileceğine ilişkin kodu veya ayrıntıları göstermezler.

Zaten verilere sahipseniz, süreç oldukça kolaydır:

R'de daewr paketi, modeli mevcut verilere uydurmanın bir örneği olarak kullanmak için güzel bir veri kümesine sahiptir.

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 

Şimdi varyansı ayarlamak ve gözlemleri simüle etmek istiyorum (daha sonra yukarıdaki analizi çalıştırın ve girdilerle karşılaştırın). Sorum şu, eğer ilgilendiğim tek şey varyansları ayarlamaksa, modeli gözlem oluşturmak için nasıl kullanabilirim? Referans makalesinde, varyans sigma ^ 2: N (0, sigma ^ 2) ile tüm rastgele etkilerin sıfır olduğunu varsayarlar. Bunun sadece rnorm (60, 0, var ^ .5) yapmak ve ardından etkileşim terimi nedeniyle terimleri eklemek kadar basit olduğunu düşünmüyorum. Etkileşim terimi kafamı karıştırıyor. Analizi çalıştırdığımda gerçek varyans bileşenlerinin makul bir tahminini alabilmem için etkileşimin rastgele etkilerle uygun şekilde hizalandığından emin olmak için bir grup matris matematiğine ihtiyacım var mı? Yoksa bundan daha mı basit?

Sağlayabileceğiniz her türlü yardım için teşekkür ederiz.

Answer 1

Temelde doğru yoldasınız.

Bunun sadece rnorm (60, 0, var ^ .5) yapmak ve ardından etkileşim terimi nedeniyle terimleri eklemek kadar basit olduğunu düşünmüyorum.

Doğru, bu yüzden sadece etkileşimin varyansı için simülasyon yapmalısın.

Karma bir model için verileri simüle etmenin en kolay yolunun model matrisini kullanmak olduğunu görüyorum. $Z$rastgele etki için. Karışık bir model için genel denklemin şu olduğunu unutmayın:

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

Ancak burada sabit bir etkimiz yok, bu yüzden sadece:

$$ Y = Zb+e $$

nerede $Z$ model matrisi rastgele etkiler ve $b$ rastgele etkiler katsayıları vektörüdür

Sorun şu ki, rastgele yapı çok basit olmadıkça, oluşturmak oldukça sıkıcı olabilir. $Z$elle. Ancak ne mutlu ki, kolay bir çözüm var - bırakın bunu sizin için yazılım yapsın. İşte sorunuzdaki model çıktısına karşılık gelen verileri kullanan bir örnek.

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

Yani burada faktörler için veri çerçevesini oluşturduk ve bir Y değişkeni oluşturmak için ona tamamen rastgele gürültü ekledik lFormulave lme4paketten formülü modele uydurmaya çalışmadan verilere karşı işlemek için kullandık . Bu işlem sırasında $ Z $ model matrisi oluşturulur ve bunun tersi $ Zt $ sonuç nesnesinde saklanır, bu nedenle oradaki son satır $ Z $ elde etmek için onu transpoze eder .

Şimdi, rastgele kesişimler için 4, 3 ve 2'nin standart sapmalarını kullandığım rastgele efektlerin kendilerini simüle ediyoruz.

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)  

Bunların girmesi gereken sırayı kontrol etmeliydim. Dokümantasyonda bunun için bazı kurallar var, ancak kodu 2 operve 2 partile çalıştırdım ve tam bir lmermodel çalıştırdım, sonra rastgele efektleri çıkarıp ranef()karşılaştırdım. getME(mymodel, "b"). Bu kafa karıştırıcıysa bana bildirin ve ben de bunun için kodu ve çıktıyı ekleyeceğim.

Ardından sonucu simüle ederiz (1 birim düzeyinde varyansla) ve lmermodeli uydururuz :

> dt$Y <- 10 + Z %*% b + rnorm(nrow(dt))
> lmer(eval(myformula), data = dt ) %>% summary()
Linear mixed model fit by REML ['lmerMod']
Formula: Y ~ (1 | part) + (1 | oper) + (1 | part:oper)
   Data: dt

REML criterion at convergence: 3776.8

Scaled residuals: 
     Min       1Q   Median       3Q      Max 
-2.42747 -0.46098  0.01696  0.46941  2.44928 

Random effects:
 Groups    Name        Variance Std.Dev.
 part:oper (Intercept) 16.833   4.103   
 oper      (Intercept) 10.183   3.191   
 part      (Intercept)  4.840   2.200   
 Residual               1.009   1.005   

Ve 4, 3, 2 ve 1 parametrelerini varyans bileşenleri olarak güzel bir şekilde kurtardığımızı görüyoruz.