Моделирование наблюдений для двухфакторного дисперсионного анализа (в 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 $ сохраняется в результирующем объекте, поэтому последняя строка там просто транспонирует ее, чтобы получить $ 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")которые сделали это очевидным. . Если это сбивает с толку, дайте мне знать, и я тоже добавлю код и вывод для этого.

Затем мы просто моделируем результат (с дисперсией на уровне единицы) и подгоняем lmerмодель: