혼합 효과 모델이있는 2- 원 분산 분석 (R)에 대한 관측치 시뮬레이션 및 실제 분산 매개 변수 복구 [Gage R & R]
R에서 Gage R & R 실험의 시뮬레이션을 만들고 싶습니다. 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
이제 분산을 설정하고 관찰을 시뮬레이션하고 싶습니다 (그런 다음 위의 분석을 실행하고 입력과 비교). 제 질문은 분산을 설정하는 데 관심이 있다면 어떻게 모델을 사용하여 관측치를 생성 할 수 있습니까? 참조 기사에서 그들은 분산 sigma ^ 2 : N (0, sigma ^ 2)를 사용하여 모든 랜덤 효과가 0이라고 가정합니다. 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 $ 가 결과 객체에 저장되므로 마지막 줄은 $ 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을 훌륭하게 복구했습니다.