混合効果モデルを使用した2元配置分散分析（R）の観測値のシミュレーションと、真の分散パラメーターの回復[ゲージR＆R]
RでのゲージR＆R実験のシミュレーションを作成したいと思います。ゲージR＆Rは、全体的な分散に対するいくつかの要因の分散の寄与を分析するために設計された実験です。多くの場合、コンテキストは、測定システムがオペレーター間の変動、部品間の変動、およびランダムな変動（再現性）の変動に起因する変動の量を知りたい測定システムです。このタイプの実験からの観測値は、通常、部分のランダム効果、オペレーターの1つ、パーツとオペレーターの交互作用、およびランダム誤差項を持つ混合効果モデルを使用してモデル化されます。各オペレーターが同じ部品を繰り返し測定することに注意してください。
私はシミュレーションを説明複製しようとしているHERE私たちは、各要因の分散を指定する場合には、観測結果を生成し、モデルに適合し、分散成分の推定値が真と比較する方法を参照してください。これらは一般的なプロセスを示していますが、差異が指定された後のデータの生成方法のコードや詳細は示していません。
すでにデータがある場合、プロセスは非常に簡単です。
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）であると想定しています。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)
これらが入る順序を確認する必要がありました。ドキュメントにはいくつかのルールがありますが、2operと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を適切に復元したことがわかります。