三层嵌套混合效应模型

Question

三层嵌套混合效应模型

6

我将尝试在rjags中建模一个三级嵌套线性混合效应模型（三级：多个个体在多个组内有多个观察结果）。在这些组中存在独特的个体集合。

相当的lme4模型为：

lmer(yN ~ x + (1 |group/indiv), data=qq)

或者

lmer(yN ~ x + (1 |group) + (1|indiv), data=qq)

我的问题是：请问如何在rjags中编写这个模型。

以下是我尝试使用rjags编写的代码，它可以编译和执行，但是个体水平的随机效应似乎被惩罚得太多了——足以表明它编码错误。

st <- "
model {

  for(i in 1:n){
      mu[i] <- beta[1] + b1[ind[i]] + b2[group[i]] + beta[2]* x[i] 
      y[i] ~ dnorm(mu[i], tau)
  }

  for(i in 1:2){  beta[i] ~ dnorm(0, 0.0001)  }

  tau ~ dgamma(0.01, 0.01)
  sigma <- sqrt(1/tau) 

  # hierarchical model
  for (i in 1:nInd) { b1[i] ~ dnorm(0, tau0) }
  for (i in 1:nGrp) { b2[i] ~ dnorm(0, tau1) }

  tau0 ~ dgamma(0.001, 0.001)
  sigma0 <- sqrt(1/tau0) 
  tau1 ~ dgamma(0.001, 0.001)
  sigma1 <- sqrt(1/tau1) 
}
"

运行模型

library(rjags)

mod <- jags.model( textConnection(st),
                 data=list(y=qq$yN, 
                           x=qq$x, 
                           ind=qq$indiv, 
                           group=qq$group,
                           n=nrow(qq),
                           nInd=length(unique(qq$indiv)),
                           nGrp=length(unique(qq$group))),
                 n.adapt=1e6,
                 inits=list(.RNG.seed=1,
                            .RNG.name="base::Wichmann-Hill")
                )
mod <- coda.samples(mod, 
                   variable.names=c("beta","b1", "b2", "sigma", "sigma0", "sigma1"), 
                   n.iter=1e6, 
                   thin=5)

summary(mod)

qq <- structure(list(yN = c(3.51, 5.13, 5.2, 7.46, 5.64, 5.14, 6.84, 
7.19, 7.77, 6, 10.97, 9.75, 5.43, 1.11, 10.31, 5.3, 4.52, 4.62, 
3.97, 4.31, 8.2, 7.24, 6.75, 0, 7.77, 4.25, 5.29, 2.46, 4.3, 
6.67, 8.72, 7.52, 6.12, 6.02, 1.48, 4.65, 7.52, 5.88, 6.06, 5.27, 
6.04, 5.36, 7.34, 6.39, 2.84, 3.95, 8.07, 7.22, 4.78, 9.92, 5.85, 
2.75, 6.34, 2.62, 7.3, 15.45, 5, 1.52, 8.3, 6.25, 16.32, 5.67, 
8.55, 5.72, 2.8, 6.06, 1.3, 11.74, 7.02, 12.85, 6.46, 3.68, 8.48, 
0.28, 0.92), x = c(-0.63, 0.18, -0.84, 1.6, 0.33, -0.82, 0.49, 
0.74, 0.58, -0.31, 1.51, 0.39, -0.62, -2.21, 1.12, -0.04, -0.02, 
0.94, 0.82, 0.59, 0.92, 0.78, 0.07, -1.99, 0.62, -0.06, -0.16, 
-1.47, -0.48, 0.42, 1.36, -0.1, 0.39, -0.05, -1.38, -0.41, -0.39, 
-0.06, 1.1, 0.76, -0.16, -0.25, 0.7, 0.56, -0.69, -0.71, 0.36, 
0.77, -0.11, 0.88, 0.4, -0.61, 0.34, -1.13, 1.43, 1.98, -0.37, 
-1.04, 0.57, -0.14, 2.4, -0.04, 0.69, 0.03, -0.74, 0.19, -1.8, 
1.47, 0.15, 2.17, 0.48, -0.71, 0.61, -0.93, -1.25), indiv = structure(c(1L, 
1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 3L, 3L, 3L, 3L, 3L, 4L, 4L, 
4L, 4L, 4L, 5L, 5L, 5L, 5L, 5L, 6L, 6L, 6L, 6L, 6L, 7L, 7L, 7L, 
7L, 7L, 8L, 8L, 8L, 8L, 8L, 9L, 9L, 9L, 9L, 9L, 10L, 10L, 10L, 
10L, 10L, 11L, 11L, 11L, 11L, 11L, 12L, 12L, 12L, 12L, 12L, 13L, 
13L, 13L, 13L, 13L, 14L, 14L, 14L, 14L, 14L, 15L, 15L, 15L, 15L, 
15L), .Label = c("a", "b", "c", "d", "e", "f", "g", "h", "i", 
"j", "k", "l", "m", "n", "o"), class = "factor"), group = structure(c(1L, 
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 
2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 3L, 3L, 3L, 
3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 4L, 4L, 4L, 4L, 
4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 5L, 5L, 5L, 5L, 5L, 
5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L), .Label = c("A", "B", 
"C", "D", "E"), class = "factor")), .Names = c("yN", "x", "indiv", 
"group"), row.names = c(NA, -75L), class = "data.frame")

在一个类似的例子中，可以通过创建交互变量并将其用作分组变量来考虑数据的嵌套结构（与先前的唯一集合组内相似）。

data(Pastes, package="lme4")

lmer(strength ~ 1 + (1|batch/cask), data=Pastes)
lmer(strength ~ 1 + (1|batch) + (1|batch:cask), data=Pastes) # equivalent

如何在 jags 中编写此代码，是否可以在不创建中间交互变量的情况下完成？

- user2957945

这只是因为JAGS和lme4的报告方式不同吗？JAGS将这些系数报告为与截距的估计差异，而lme4为该随机效应提供了估计的截距值。如果您将beta [1]添加回到每个b1中，则估计值应该非常接近lme4中报告的值。 - mfidino

@M_Fidino；不，我认为不是这样。我已经设置了JAGS代码，以复制lme4呈现结果的方式（包括固定效应和随机偏差）。b预测应该与ranef（lme4model）相同。 - user2957945

@M_Fidino；也许仅仅是先验选择的症状——我不确定我是否正确地编写了语法以解决嵌套结构问题。（感谢您的查看） - user2957945

JAGS模型假设随机效应之间没有协方差（因为它们各自从自己的正态分布中抽取），而lme4模型则有。也许这就是造成差异的原因？ - mfidino

@M_Fidino；也许，谢谢，然而，如果我在开头运行lmer模型中的任何一个，定义或返回的随机项之间没有相关性（如果我理解正确）。 - user2957945

1个回答

网页内容由stack overflow 提供, 点击上面的

可以查看英文原文，
原文链接

- mfidino · Accepted Answer

对于嵌套效应，您需要将每个效应与其所在的特定组链接起来。当前的 JAGS 模型目前不支持此功能。为此，您需要另一个向量，将个体与组链接起来。

unq_ind_group <- qq[,3:4]
unq_ind_group <- unq_ind_group[!duplicated(unq_ind_group),]

更新后的模型：

st <- "
model {
for(i in 1:n){
mu[i] <- beta[1] + b1[ind[i]] + b2[group[i]] + beta[2]* x[i] 
y[i] ~ dnorm(mu[i], tau)
}
for(i in 1:2){  beta[i] ~ dnorm(0, 0.0001)  }
tau ~ dgamma(0.01, 0.01)
sigma <- sqrt(1/tau) 
# hierarchical model
for (i in 1:nGrp) { b2[i] ~ dnorm(0, tau1) }
for (i in 1:nInd) { b1[i] ~ dnorm(b2[ind_per_group[i]], tau0) }
tau0 ~ dgamma(0.001, 0.001)
sigma0 <- sqrt(1/tau0) 
tau1 ~ dgamma(0.001, 0.001)
sigma1 <- sqrt(1/tau1) 
}
"
# fit the model
mod <- jags.model( textConnection(st),
        data=list(y=qq$yN, 
        x=qq$x, 
        ind=qq$indiv, 
        group=qq$group,
        ind_per_group = unq_ind_group$group,
        n=nrow(qq),
        nInd=length(unique(qq$indiv)),
        nGrp=length(unique(qq$group))),
        n.adapt=1e6,
        inits=list(.RNG.seed=1,
       .RNG.name="base::Wichmann-Hill")
)

mod <- coda.samples(mod, 
    variable.names=c("beta","b1", "b2", "sigma", "sigma0", "sigma1"), 
    n.iter=1e6, 
    thin=5)

这里是上述模型和lme4中嵌套模型标准差的比较。

m2 <- lmer(yN ~ x + (1 |group/indiv), data=qq)
summary(m2)

这个模型总结告诉我们：

单体：组的标准偏差为0.7909
组标准偏差为0
残差偏差为1.3629

这是一个比较不同模型间估计值的图表。白色点是JAGS估计值，黑色点是lme4估计值，竖线是JAGS的95%可信区间。

此外，你设置的随机效应精度项的先验分布大部分集中在零点，这会影响后验分布。这是因为每组之间只有少数个体，所以数据无法超过先验分布。注意，sigma0的可信区间是三者中最大的，反映了这一估计的不确定性。如果你的目标是更接近lme4的估计值，可以设置dgamma(0.1,0.1)先验。

更新：

这是一个将JAGS模型的随机效应与lme4进行比较的图表。与前面的图表类似，白色点是JAGS的中位数估计值，黑色点是通过ranef(m2)从lme4获取的估计值，竖线是JAGS的95%可信区间。从这个图表可以看出，由于sigma0被估计为较小的值，所有随机效应的JAGS估计值都被拉向零点。

这是我如何修改JAGS模型以跟踪这些随机效应作为衍生参数的方式。然后，我只需将"b_pred"作为coda.samples中要跟踪的另一个元素添加到variable.names参数中。

st <- "
model {

for(i in 1:n){
mu[i] <- beta[1] + b1[ind[i]] + b2[group[i]] + beta[2]* x[i] 
y[i] ~ dnorm(mu[i], tau)
}

for(i in 1:2){  beta[i] ~ dnorm(0, 0.0001)  }

tau ~ dgamma(0.01, 0.01)
sigma <- sqrt(1/tau) 

# hierarchical model
for (i in 1:nGrp) { b2[i] ~ dnorm(0, tau1) }
for (i in 1:nInd) { b1[i] ~ dnorm(b2[ind_per_group[i]], tau0) }

tau0 ~ dgamma(0.001, 0.001)
sigma0 <- sqrt(1/tau0) 
tau1 ~ dgamma(0.001, 0.001)
sigma1 <- sqrt(1/tau1) 
# calculate random effects
for(i in 1:nInd) {b_pred[i] <- b1[i] + b2[ind_per_group[i]]}

}
"