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() 可以使随机效应相关（参见第 3.4 节和 Hadfield's 2010 jstatsoft-article 中的 table 2）

首先，由于您只有一个因变量 (y)，因此先验中的 G 部分（参见等式 4 和第 3.6 节 Hadfield's 2010 jstatsoft-article) 对于随机效应方差，只需要有一个名为 G1 的列表元素。此列表元素不是实际的先验分布 - Hadfield 将其指定为 inverse-Wishart 分布。但是使用 G1 您可以指定此逆惠沙特分布的参数，即比例矩阵（维基百科符号中的 $\Psi$ 和 MCMCglmm 符号中的 V ）和自由度（维基百科符号中的 $\nu$ 和 MCMCglmm 符号中的 nu）。由于您有两个随机效应（截距和斜率）V 必须是 2 x 2 矩阵。一个常见的选择是二维单位矩阵 diag(2)。哈德菲尔德经常使用 nu = 0.002 作为自由度（参见他的课程笔记）

现在，您还必须在先验中为残差方差指定 R 部分。此处再次由 Hadfield 指定了逆 Whishart 分布，让用户指定其参数。因为我们只有一个残差方差，所以 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)

在 mod0$Sol 中找到固定效应的 Gibbs 采样器抽取，在 mod0$VCV.

中找到方差参数的抽取

通常二项式模型要求残差方差固定，所以我们设置残差固定在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)

通过比较 mod0$VCV[, 5] 和 mod1$VCV[, 5] 可以看出差异。在后一种情况下，所有条目都是 0.5 指定的。

MCMCglmm 二项式模型先验

MCMCglmm binomial model prior

r

mcmc

glm

mixed-models