三级嵌套混合效应模型
Three-level nested mixed-effects model
我正在尝试在 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
中编码,是否可以在不创建中间交互变量的情况下完成?
对于嵌套效果,您需要 link 将单个效果分配给它们所在的特定组。当前的 JAGS 模型当前不执行此操作。为此,您需要另一个向量,该向量 link 将个体与群体联系起来。
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)
这个模型的总结告诉我们
- 与 indiv:group 的标准差是 0.7909
- 组标准差为0
- 残差为1.3629
这是一个比较模型之间估计值的图。白点是 JAGS 估计值,黑点来自 lme4,垂直线是 JAGS 的 95% 可信区间。
此外,您为随机效应的精度项设置的先验大部分质量为零,这将影响后验分布。这是因为每个组中的个体很少,因此数据不会超过先前的数据。请注意,sigma0
的可信区间是三个区间中最大的,这反映了该估计的不确定性。设置更接近 lme4 的 dgamma(0.1,0.1)
先验 returns 估计值(如果这是您的目标)。
更新:
这是比较 JAGS
模型和 lme4
模型的随机效应的图表。和之前的剧情一样。白点是 JAGS
的中值估计值,黑点是 lme4
通过 ranef(m2)
的估计值,垂直线是 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]]}
}
"
我正在尝试在 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
中编码,是否可以在不创建中间交互变量的情况下完成?
对于嵌套效果,您需要 link 将单个效果分配给它们所在的特定组。当前的 JAGS 模型当前不执行此操作。为此,您需要另一个向量,该向量 link 将个体与群体联系起来。
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)
这个模型的总结告诉我们
- 与 indiv:group 的标准差是 0.7909
- 组标准差为0
- 残差为1.3629
这是一个比较模型之间估计值的图。白点是 JAGS 估计值,黑点来自 lme4,垂直线是 JAGS 的 95% 可信区间。
此外,您为随机效应的精度项设置的先验大部分质量为零,这将影响后验分布。这是因为每个组中的个体很少,因此数据不会超过先前的数据。请注意,sigma0
的可信区间是三个区间中最大的,这反映了该估计的不确定性。设置更接近 lme4 的 dgamma(0.1,0.1)
先验 returns 估计值(如果这是您的目标)。
更新:
这是比较 JAGS
模型和 lme4
模型的随机效应的图表。和之前的剧情一样。白点是 JAGS
的中值估计值,黑点是 lme4
通过 ranef(m2)
的估计值,垂直线是 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]]}
}
"