2-तरफ़ा एनोवा (R) w / मिश्रित प्रभाव मॉडल के लिए टिप्पणियों का अनुकरण करना और वास्तविक विचरण मापदंडों को ठीक करना [Gage R & R]
मैं R. A Gage R & R में एक Gage R & R प्रयोग का अनुकरण बनाना चाहूंगा, जो एक ऐसा डिज़ाइन है जिसे समग्र भिन्नता के सापेक्ष कई कारकों के भिन्नता योगदान का विश्लेषण करने के लिए डिज़ाइन किया गया है। संदर्भ अक्सर माप प्रणाली है जहां हम यह जानना चाहते हैं कि ऑपरेटर-से-ऑपरेटर भिन्नता, भाग-से-भाग भिन्नता, और यादृच्छिक भिन्नता (repeatability) भिन्नता के कारण माप प्रणाली कितनी भिन्नता है। इस तरह के प्रयोगों से टिप्पणियों को आमतौर पर मिश्रित प्रभाव वाले मॉडल का उपयोग करके भाग के लिए यादृच्छिक प्रभाव, ऑपरेटर के लिए एक, एक भाग: ऑपरेटर इंटरैक्शन और एक यादृच्छिक त्रुटि शब्द के साथ तैयार किया जाता है। ध्यान दें कि प्रत्येक ऑपरेटर एक ही हिस्से की बार-बार माप करता है।
मैं यहां वर्णित सिमुलेशन को दोहराने की कोशिश कर रहा हूं, जहां हम प्रत्येक कारक के लिए विचरण निर्दिष्ट करते हैं, अवलोकनों को उत्पन्न करते हैं, फिर एक मॉडल फिट करते हैं और देखते हैं कि कैसे विचरण घटकों के अनुमान सच के साथ तुलना करते हैं। वे सामान्य प्रक्रिया दिखाते हैं लेकिन कोड या बारीकियों को निर्दिष्ट करने के लिए नहीं कि डेटा कैसे उत्पन्न करें।
यदि आपके पास पहले से ही डेटा है, तो प्रक्रिया बहुत आसान है:
आर में, डेवेर पैकेज में मौजूदा डेटा के लिए मॉडल को फिट करने के उदाहरण के रूप में उपयोग करने के लिए एक अच्छा डेटासेट है
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, sigma ^ 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
इसलिए यहां हमने सिर्फ कारकों के लिए डेटा फ्रेम बनाया है, और वाई चर बनाने के लिए इसमें विशुद्ध रूप से यादृच्छिक शोर जोड़ा है और मॉडल को फिट करने का प्रयास किए बिना डेटा के खिलाफ सूत्र को संसाधित करने के lFormulaलिए lme4पैकेज से उपयोग किया जाता है । इस प्रक्रिया के दौरान $ Z $ मॉडल मैट्रिक्स का निर्माण किया जाता है और इसका उलटा $ Zt $ परिणामी वस्तु में जमा हो जाता है, इसलिए अंतिम पंक्ति बस इसे $ 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 के रूप में पैरामीटर 4, 3, 2 और 1 को पुनर्प्राप्त किया है