使用 JAGS for R 的分层模型显示错误的均值数
Hierarchical model using JAGS for R shows wrong number of means
我有 2 列数据,y
和 grp
,我正在尝试创建如上所示的 JAGS
模型。 grp
是组,我有 5 个组。以下代码来自here。我使用这段代码是因为标题 Model and Data
下的描述看起来几乎像这个分层模型。
但我在查看摘要时只得到一个 mu
。应该有 5 mu's
,每组一个。有人可以更正代码吗?您也可以指出其他地方提供的类似示例,我可以尝试对其进行修改。我在代码中遗漏了一些东西,我相信代码可能类似于 question 但是当我这样修改它时,即使有 5 种方法,我似乎也没有得到正确的方法。
不确定这个问题是否属于 math stackexchange。
mod_string = " model {
for (i in 1:length(y) {
theta[i] ~ dnorm(mu[grp[i]], invTau2)
y[i] ~ dnorm(theta[i], 1/sig)
}
mu ~ dnorm(0, 1e6)
invTau2 ~ dgamma(1.0/2.0, 1.3/2.0)
tau2 <- 1/invTau2
invgamma2 ~ dgamma(1.0/2.0, 2.1/2.0)
sig = 1/invgamma2
} "
summary(mod_sim)
Iterations = 2001:52000
Thinning interval = 1
Number of chains = 3
Sample size per chain = 50000
1. Empirical mean and standard deviation for each variable,
plus standard error of the mean:
Mean SD Naive SE Time-series SE
mu 5.639e-07 0.001 2.582e-06 2.582e-06
sig 1.570e+00 1.888 4.874e-03 7.068e-03
您似乎在使用嵌套索引,其中 grp
是一个长度为 y
的向量,描述了每个元素 y
属于哪个组。如果每个组都需要它自己的平均值,那么您需要生成 5 theta
个先验。这在它自己的循环中最容易完成。因此,您的模型字符串看起来像:
mod_string <- " model {
# loop for y
for (i in 1:length(y)) {
y[i] ~ dnorm(theta[grp[i]], sig)
}
# loop for theta
for (g in 1:ntheta){
theta[g] ~ dnorm(mu, tau2)
}
# other priors. JAGS uses precision so you need to use the reciprocal
# of 1e6 for the variance of mu.
mu ~ dnorm(0, 0.000001)
invTau2 ~ dgamma(1.0/2.0, 1.3/2.0)
tau2 <- 1/invTau2
invgamma2 ~ dgamma(1.0/2.0, 2.1/2.0)
sig <- 1/invgamma2
} "
请注意,您还需要向模型提供标量 ntheta
,这将是唯一组的数量(在本例中为 ntheta = 5
)。现在,如果您跟踪 theta
,您将得到 5 个估计值,而 mu
有 1 个估计值(这是所有组的总估计平均值)。
我有 2 列数据,y
和 grp
,我正在尝试创建如上所示的 JAGS
模型。 grp
是组,我有 5 个组。以下代码来自here。我使用这段代码是因为标题 Model and Data
下的描述看起来几乎像这个分层模型。
但我在查看摘要时只得到一个 mu
。应该有 5 mu's
,每组一个。有人可以更正代码吗?您也可以指出其他地方提供的类似示例,我可以尝试对其进行修改。我在代码中遗漏了一些东西,我相信代码可能类似于 question 但是当我这样修改它时,即使有 5 种方法,我似乎也没有得到正确的方法。
不确定这个问题是否属于 math stackexchange。
mod_string = " model {
for (i in 1:length(y) {
theta[i] ~ dnorm(mu[grp[i]], invTau2)
y[i] ~ dnorm(theta[i], 1/sig)
}
mu ~ dnorm(0, 1e6)
invTau2 ~ dgamma(1.0/2.0, 1.3/2.0)
tau2 <- 1/invTau2
invgamma2 ~ dgamma(1.0/2.0, 2.1/2.0)
sig = 1/invgamma2
} "
summary(mod_sim)
Iterations = 2001:52000
Thinning interval = 1
Number of chains = 3
Sample size per chain = 50000
1. Empirical mean and standard deviation for each variable,
plus standard error of the mean:
Mean SD Naive SE Time-series SE
mu 5.639e-07 0.001 2.582e-06 2.582e-06
sig 1.570e+00 1.888 4.874e-03 7.068e-03
您似乎在使用嵌套索引,其中 grp
是一个长度为 y
的向量,描述了每个元素 y
属于哪个组。如果每个组都需要它自己的平均值,那么您需要生成 5 theta
个先验。这在它自己的循环中最容易完成。因此,您的模型字符串看起来像:
mod_string <- " model {
# loop for y
for (i in 1:length(y)) {
y[i] ~ dnorm(theta[grp[i]], sig)
}
# loop for theta
for (g in 1:ntheta){
theta[g] ~ dnorm(mu, tau2)
}
# other priors. JAGS uses precision so you need to use the reciprocal
# of 1e6 for the variance of mu.
mu ~ dnorm(0, 0.000001)
invTau2 ~ dgamma(1.0/2.0, 1.3/2.0)
tau2 <- 1/invTau2
invgamma2 ~ dgamma(1.0/2.0, 2.1/2.0)
sig <- 1/invgamma2
} "
请注意,您还需要向模型提供标量 ntheta
,这将是唯一组的数量(在本例中为 ntheta = 5
)。现在,如果您跟踪 theta
,您将得到 5 个估计值,而 mu
有 1 个估计值(这是所有组的总估计平均值)。