การจำลองการสังเกตสำหรับ 2-way ANOVA (ใน R) ที่มีโมเดลเอฟเฟกต์ผสมและการกู้คืนพารามิเตอร์ความแปรปรวนที่แท้จริง [Gage R&R]

Aug 18 2020

ฉันต้องการสร้างแบบจำลองของการทดลอง Gage R&R ใน R A Gage R&R เป็นการทดลองที่ออกแบบมาเพื่อวิเคราะห์การสนับสนุนความแปรปรวนของปัจจัยหลายอย่างที่สัมพันธ์กับความแปรปรวนโดยรวม บริบทมักเป็นระบบการวัดที่เราต้องการทราบว่าระบบการวัดมีความผันแปรมากน้อยเพียงใดเนื่องจากรูปแบบของตัวดำเนินการต่อตัวดำเนินการรูปแบบส่วนต่อส่วนและการเปลี่ยนแปลงแบบสุ่ม (ความสามารถในการทำซ้ำ) โดยทั่วไปการสังเกตจากการทดลองประเภทนี้จะสร้างแบบจำลองโดยใช้แบบจำลองเอฟเฟกต์ผสมที่มีเอฟเฟกต์แบบสุ่มสำหรับส่วนหนึ่งสำหรับตัวดำเนินการส่วนหนึ่ง: การโต้ตอบของตัวดำเนินการและเงื่อนไขข้อผิดพลาดแบบสุ่ม โปรดทราบว่าตัวดำเนินการแต่ละตัวทำการวัดซ้ำในส่วนเดียวกัน

ฉันกำลังพยายามจำลองแบบจำลองที่อธิบายไว้ที่นี่ซึ่งเราระบุความแปรปรวนของแต่ละปัจจัยสร้างการสังเกตจากนั้นให้พอดีกับแบบจำลองและดูว่าค่าประมาณของส่วนประกอบความแปรปรวนเปรียบเทียบกับความจริงอย่างไร แสดงกระบวนการทั่วไป แต่ไม่แสดงรหัสหรือเฉพาะสำหรับวิธีสร้างข้อมูลเมื่อระบุความแปรปรวน

หากคุณมีข้อมูลอยู่แล้วกระบวนการนี้ค่อนข้างง่าย:

ใน R แพ็กเกจ daewr มีชุดข้อมูลที่ดีเพื่อใช้เป็นตัวอย่างในการปรับโมเดลเข้ากับข้อมูลที่มีอยู่

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 

ตอนนี้ฉันต้องการตั้งค่าความแปรปรวนและจำลองการสังเกต (จากนั้นเรียกใช้การวิเคราะห์ด้านบนและเปรียบเทียบกับอินพุต) คำถามของฉันคือฉันจะใช้แบบจำลองเพื่อสร้างข้อสังเกตได้อย่างไรหากสิ่งที่ฉันสนใจคือการตั้งค่าความแปรปรวน ในบทความอ้างอิงพวกเขาถือว่าเอฟเฟกต์แบบสุ่มทั้งหมดเป็นศูนย์โดยมีซิกม่าความแปรปรวน ^ 2: N (0, ซิกม่า ^ 2) ฉันไม่คิดว่ามันง่ายเหมือนแค่ทำ rnorm (60, 0, var ^ .5) จากนั้นเพิ่มเงื่อนไขเนื่องจากเงื่อนไขการโต้ตอบ คำโต้ตอบทำให้ฉันสับสน ฉันต้องการคณิตศาสตร์เมทริกซ์จำนวนมากเพื่อให้แน่ใจว่าการโต้ตอบสอดคล้องกับเอฟเฟกต์แบบสุ่มอย่างเหมาะสมหรือไม่เช่นนั้นเมื่อฉันทำการวิเคราะห์ฉันจะได้ค่าประมาณที่สมเหตุสมผลของส่วนประกอบความแปรปรวนจริงหรือไม่ หรือว่าง่ายกว่านั้น?

ขอบคุณสำหรับความช่วยเหลือที่คุณสามารถให้ได้

คำตอบ

5 RobertLong Aug 19 2020 at 13:02

คุณมาถูกทางแล้ว

ฉันไม่คิดว่ามันง่ายเหมือนแค่ทำ rnorm (60, 0, var ^ .5) จากนั้นเพิ่มเงื่อนไขเนื่องจากเงื่อนไขการโต้ตอบ

ถูกต้องดังนั้นคุณต้องจำลองความแปรปรวนของการโต้ตอบด้วย

ฉันพบว่าวิธีที่ง่ายที่สุดในการจำลองข้อมูลสำหรับแบบจำลองผสมคือการใช้เมทริกซ์แบบจำลอง $Z$สำหรับเอฟเฟกต์แบบสุ่ม จำไว้ว่าสมการทั่วไปสำหรับแบบจำลองผสมคือ:

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

แต่ที่นี่เราไม่มีเอฟเฟกต์คงที่ดังนั้นจึงเป็นเพียง:

$$ Y = Zb+e $$

ที่ไหน $Z$ คือเมทริกซ์แบบจำลองเอฟเฟกต์แบบสุ่มและ $b$ คือเวกเตอร์สัมประสิทธิ์เอฟเฟกต์สุ่ม

ปัญหาคือเว้นแต่โครงสร้างแบบสุ่มจะง่ายมากการสร้างก็ค่อนข้างน่าเบื่อ $Z$ด้วยมือ. แต่อย่างมีความสุขมีวิธีง่ายๆเพียงแค่ปล่อยให้ซอฟต์แวร์ทำเพื่อคุณ นี่คือตัวอย่างการใช้ข้อมูลที่ตรงกับผลลัพธ์ของโมเดลในคำถามของคุณ

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

ดังนั้นที่นี่เราเพิ่งสร้างกรอบข้อมูลสำหรับปัจจัยและเพิ่มสัญญาณรบกวนแบบสุ่มอย่างหมดจดเพื่อสร้างตัวแปร Y และใช้lFormulaจากlme4แพ็คเกจเพื่อประมวลผลสูตรกับข้อมูลโดยไม่พยายามปรับให้พอดีกับโมเดล ในช่วงนี้การประมวลผล$ Z $รุ่นเมทริกซ์ถูกสร้างและมันผกผัน$ $ Ztจะถูกเก็บไว้ในวัตถุที่เกิดขึ้นเพื่อให้บรรทัดสุดท้ายมีเพียงแค่ transposes จะได้รับ$ Z $

ตอนนี้เราจำลองเอฟเฟกต์แบบสุ่มด้วยตัวเองโดยที่ฉันใช้ค่าเบี่ยงเบนมาตรฐาน 4, 3 และ 2 สำหรับการสกัดกั้นแบบสุ่ม

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)  

ฉันต้องตรวจสอบลำดับว่าสิ่งเหล่านี้ควรเข้าอย่างไรมีกฎบางอย่างสำหรับสิ่งนี้ในเอกสาร แต่ฉันเพียงแค่รันโค้ดด้วย 2 operและ 2 partและรันแบบlmerจำลองเต็มรูปแบบจากนั้นแยกเอฟเฟกต์แบบสุ่มด้วยranef()และเปรียบเทียบกับสิ่งgetME(mymodel, "b")ที่ทำให้ชัดเจน . หากสิ่งนี้ทำให้สับสนโปรดแจ้งให้เราทราบและฉันจะเพิ่มโค้ดและเอาต์พุตสำหรับสิ่งนั้นด้วย

จากนั้นเราเพียงแค่จำลองผลลัพธ์ (โดยมีความแปรปรวนระดับหน่วยเป็น 1) และพอดีกับlmerโมเดล:

> 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   

และเราเห็นว่าเราได้กู้คืนพารามิเตอร์ 4, 3, 2 และ 1 เป็นส่วนประกอบความแปรปรวนแล้ว