การจำลองการสังเกตสำหรับ 2-way ANOVA (ใน R) ที่มีโมเดลเอฟเฟกต์ผสมและการกู้คืนพารามิเตอร์ความแปรปรวนที่แท้จริง [Gage R&R]
ฉันต้องการสร้างแบบจำลองของการทดลอง 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) จากนั้นเพิ่มเงื่อนไขเนื่องจากเงื่อนไขการโต้ตอบ คำโต้ตอบทำให้ฉันสับสน ฉันต้องการคณิตศาสตร์เมทริกซ์จำนวนมากเพื่อให้แน่ใจว่าการโต้ตอบสอดคล้องกับเอฟเฟกต์แบบสุ่มอย่างเหมาะสมหรือไม่เช่นนั้นเมื่อฉันทำการวิเคราะห์ฉันจะได้ค่าประมาณที่สมเหตุสมผลของส่วนประกอบความแปรปรวนจริงหรือไม่ หรือว่าง่ายกว่านั้น?
ขอบคุณสำหรับความช่วยเหลือที่คุณสามารถให้ได้
คำตอบ
คุณมาถูกทางแล้ว
ฉันไม่คิดว่ามันง่ายเหมือนแค่ทำ 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 เป็นส่วนประกอบความแปรปรวนแล้ว