Simulasi pengamatan untuk ANOVA 2 arah (dalam R) dengan model efek campuran dan memulihkan parameter varian sebenarnya [Gage R&R]

Aug 18 2020

Saya ingin membuat simulasi eksperimen Gage R&R di R. Gage R&R adalah eksperimen yang dirancang untuk menganalisis kontribusi varians dari beberapa faktor relatif terhadap varian keseluruhan. Konteksnya sering kali merupakan sistem pengukuran tempat kami ingin mengetahui seberapa banyak variasi sistem pengukuran karena variasi operator-ke-operator, variasi bagian-ke-bagian, dan variasi variasi acak (pengulangan). Pengamatan dari jenis eksperimen ini biasanya dimodelkan menggunakan model efek campuran dengan efek acak untuk sebagian, satu untuk operator, satu bagian: interaksi operator, dan istilah kesalahan acak. Perhatikan bahwa setiap operator melakukan pengukuran berulang untuk bagian yang sama.

Saya mencoba untuk mereplikasi simulasi yang dijelaskan DI SINI di mana kita menentukan varians untuk setiap faktor, menghasilkan pengamatan, kemudian menyesuaikan model dan melihat bagaimana perkiraan komponen varians dibandingkan dengan yang sebenarnya. Mereka menunjukkan proses umum tetapi bukan kode atau spesifik untuk bagaimana menghasilkan data setelah varians ditentukan.

Jika sudah punya datanya, prosesnya cukup mudah:

Di R, paket daewr memiliki dataset yang bagus untuk digunakan sebagai contoh menyesuaikan model ke data yang ada

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

Sekarang saya ingin mengatur varians dan mensimulasikan pengamatan (kemudian menjalankan analisis di atas dan membandingkan dengan input). Pertanyaan saya adalah, bagaimana saya bisa menggunakan model untuk menghasilkan pengamatan jika yang saya pedulikan hanyalah mengatur varians? Dalam artikel referensi, mereka menganggap semua efek acak adalah nol dengan varians sigma ^ 2: N (0, sigma ^ 2). Saya tidak berpikir itu sesederhana hanya melakukan rnorm (60, 0, var ^ .5) dan kemudian menambahkan istilah karena istilah interaksi. Istilah interaksi membingungkan saya. Apakah saya memerlukan sekumpulan matematika matriks untuk memastikan interaksi selaras dengan efek acak sehingga ketika saya menjalankan analisis, saya bisa mendapatkan perkiraan yang masuk akal dari komponen varian sebenarnya? Atau lebih sederhana dari itu?

Terima kasih atas bantuan yang Anda berikan.

Jawaban

5 RobertLong Aug 19 2020 at 13:02

Anda pada dasarnya berada di jalur yang benar.

Saya tidak berpikir itu sesederhana hanya melakukan rnorm (60, 0, var ^ .5) dan kemudian menambahkan istilah karena istilah interaksi.

Benar, jadi Anda hanya perlu melakukan simulasi untuk varians interaksi juga.

Saya menemukan bahwa cara termudah untuk mensimulasikan data untuk model campuran melakukan ini adalah dengan menggunakan model matriks, $Z$untuk efek acak. Ingatlah bahwa persamaan umum untuk model campuran adalah:

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

Tapi disini kita tidak memiliki efek tetap jadi hanya saja:

$$ Y = Zb+e $$

dimana $Z$ adalah model matriks efek acak dan $b$ adalah vektor koefisien efek acak

Masalahnya adalah bahwa kecuali struktur acak sangat sederhana, pembuatannya bisa sangat membosankan $Z$dengan tangan. Namun, untungnya, ada solusi yang mudah - biarkan perangkat lunak melakukannya untuk Anda. Berikut adalah contoh penggunaan data yang sesuai dengan keluaran model dalam pertanyaan Anda.

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

Jadi di sini kita baru saja membuat kerangka data untuk faktor-faktor tersebut, dan menambahkan gangguan acak murni ke dalamnya untuk membuat variabel Y dan digunakan lFormuladari lme4paket untuk memproses rumus terhadap data tanpa mencoba menyesuaikan model. Selama pemrosesan ini, matriks model $ Z $ dibangun dan kebalikannya $ Zt $ disimpan dalam objek yang dihasilkan, jadi baris terakhir di sana hanya mentransposisi untuk mendapatkan $ Z $ .

Sekarang kami mensimulasikan efek acak itu sendiri di mana saya menggunakan deviasi standar 4, 3 dan 2 untuk penyadapan acak.

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)

Saya harus memeriksa urutan apakah ini harus masuk. Ada beberapa aturan untuk ini dalam dokumentasi tetapi saya hanya menjalankan kode dengan 2 operdan 2 partdan menjalankan lmermodel penuh kemudian mengekstrak efek acak dengan ranef()dan membandingkannya getME(mymodel, "b")yang membuatnya jelas . Jika ini membingungkan beri tahu saya dan saya akan menambahkan kode dan output untuk itu juga.