MCMCglmm二項モデル事前

Question

RパッケージMCMCglmmで二項モデルを推定したいと思います。モデルは、固定部分とランダム部分の両方として、切片と勾配を組み込むものとする。受け入れたpriorを指定するにはどうすればよいですか？ （注、here is a similar questionが、はるかに複雑な設定で）

データを想定し、以下の形式があります：実際に

y   x cluster 
1 0 -0.56047565  1 
2 1 -0.23017749  1 
3 0 1.55870831  1 
4 1 0.07050839  1 
5 0 0.12928774  1 
6 1 1.71506499  1 

---

を、データが

set.seed(123) 
nj <- 15 # number of individuals per cluster 
J <- 30 # number of clusters 
n <- nj * J 
x <- rnorm(n) 
y <- rbinom(n, 1, prob = 0.6) 
cluster <- factor(rep(1:nj, each = J)) 

dat <- data.frame(y = y, x = x, cluster = cluster) 

によって生成されています

Answer 1

モデルに関する質問の情報は、fixed = y ~ 1 + xとrandom = ~ us(1 + x):clusterと指定することを推奨します。 us()を使用すると、あなたが唯一の従属変数（y）を持っているようなランダム効果は、（参照：セクション3.4およびHadfield's 2010 jstatsoft-articleの表2）すべての

まず相関することが以前（CFでG一部を許可します（式4およびセクション3.6のHadfield's 2010 jstatsoft-article）には、G1という1つのリスト要素が必要です。このリスト要素は実際の事前分布ではなく、Hadfieldによってinverse-Wishart distributionと指定されています。しかしG1で、あなたはスケール行列（Wikipediaの表記とMCMCglmm表記でVで）と自由度（Wikipediaの表記とMCMCglmm表記でnuで）あるこの逆Whishart分布のパラメータを指定します。 2つのランダムな効果（切片と勾配）があるので、Vは2 x 2の行列でなければなりません。頻繁な選択は、二次元の単位行列diag(2)です。ハドフィールドは、多くの場合、自由度のために（参照his course notes）今

をnu = 0.002を使用して、あなたはまた、残差分散のために事前にR一部を指定する必要があります。この場合も、Hadfieldによって逆Whishart分布が指定され、ユーザはそのパラメータを指定する必要があります。残差分散が1つしかないので、Vはスカラーでなければなりません（V = 0.5と言うことができます）。 Rのオプション要素はfixです。この要素を使用して、残差分散を特定の値に固定するか（fix = TRUEまたはfix = 1を書かなければならないか）（fix = FALSEまたはfix = 0）を指定します。残差分散を0.5に固定しないことに注意してください。fix = 0.5！ですから、Hadfieldのコースノートfix = 1で見つけた場合はfix = TRUEと読んで、Vの値が修正されていることを確認してください。これにより

prior0 <- list(G = list(G1 = list(V = diag(2), nu = 0.002)), 
       R = list(V = 0.5, nu = 0.002, fix = FALSE)) 

を我々はMCMCglmmを実行することができます前：

library("MCMCglmm") # for MCMCglmm()  
set.seed(123) 

mod0 <- MCMCglmm(fixed = y ~ 1 + x, 
      random = ~ us(1 + x):cluster, 
      data = dat, 
      family = "categorical", 
      prior = prior0) 

固定効果があるためにザ・はギブスサンプラーから描く、我々は次の前のように設定すべてtogehter

mod0$Solにある場合、分散パラメータの値はmod0$VCVになります。

は、通常、二項モデルは、固定される残差分散を必要とするので、我々は、残留分散がmod1$VCV[, 5]にmod0$VCV[, 5]を比較することによって見ることができる0.5

set.seed(123) 

prior1 <- list(G = list(G1 = list(V = diag(2), nu = 0.002)), 
      R = list(V = 0.5, nu = 0.002, fix = TRUE)) 

mod1 <- MCMCglmm(fixed = y ~ 1 + x, 
      random = ~ us(1 + x):cluster, 
      data = dat, 
      family = "categorical", 
      prior = prior1) 

差に固定されるように設定しました。後者の場合、すべてのエントリは指定された通り0.5です。