r-rjags方差分析，每组具有单独的方差

Question

我正试图使用rjags对三组中的每一组进行方差分析。我正在使用"；PlantGrowth"；R.中内置的数据

这是型号：

mod_string = " model {
for (i in 1:length(y)) {
y[i] ~ dnorm(mu[grp[i]], sig[grp[i]])
}
for (j in 1:3) {
mu[j] ~ dnorm(0.0, 1.0/1.0e6)
}
for (j in 1:3) {
sig[j] ~ dnorm(0.0, 1.0/1.0e6)
}
} "
set.seed(82)
str(PlantGrowth)
data_jags = list(y=PlantGrowth$weight, 
grp=as.numeric(PlantGrowth$group))
params = c("mu", "sig")
inits = function() {
inits = list("mu"=rnorm(3,0.0,100.0), "sig"=rnorm(3,0.0,1.0))
}

此模型中的inits函数不起作用。有人能帮忙吗？

Answer 1

我建议下一种方法。首先，需要对贝叶斯模型进行训练。所以你的步伐很好。考虑到你在系数上定义了极低的先验。因此，当您将初始值设置得过高时，可能会出现与链中的收敛相关的问题。模型链很好。我只是稍微改进了初始值函数。这是代码。为了完成方差分析，你必须处理你感兴趣的参数的所有结果。这意味着计算后验可信区间等。这样，你的主要输入将是来自链的参数。这里是获得您想要的参数的代码：

library(rjags)
set.seed(82)
#Model string
mod_string = " model {
for (i in 1:length(y)) {
y[i] ~ dnorm(mu[grp[i]], sig[grp[i]])
}
for (j in 1:3) {
mu[j] ~ dnorm(0.0, 1.0/1.0e6)
}
for (j in 1:3) {
sig[j] ~ dnorm(0.0, 1.0/1.0e6)
}
} "
#Data
data(PlantGrowth)
#Jags inputs
data_jags = list(y=PlantGrowth$weight, 
grp=as.numeric(PlantGrowth$group))
params = c("mu", "sig")
#Initial values
inits = function() {
inits = list("mu"=runif(3,0,1), "sig"=runif(3,0,1))
}
#Train the model
#jags model
results.model <- jags.model(file=textConnection(mod_string),
data=data_jags, inits=inits(),
n.chains=3,quiet = T)
#Update
update(results.model, n.iter=10000,progress.bar = 'none')
#Sampling
results.model.1 <- coda.samples(results.model,
variable.names=params,n.iter=10000,
progress.bar = 'none')
#Extract a sample
d1 <- as.data.frame(do.call(rbind,results.model.1))

输出将存储在d1上。这里有几行：

head(d1)
mu[1]    mu[2]    mu[3]   sig[1]   sig[2]   sig[3]
1 5.121681 4.725264 5.601421 3.625515 3.096606 9.499369
2 5.399050 4.746094 5.551337 1.260319 3.055592 7.548953
3 4.904069 4.816449 5.655899 1.414513 2.555520 4.303972
4 4.538518 4.730499 5.481598 1.297098 1.597760 5.643879
5 4.908547 4.339875 5.452959 3.710132 1.062263 5.865650
6 5.250289 4.845275 5.554499 2.238207 2.382826 6.492258

有了这些值，你必须计算后验均值和后验可信区间，以便与方差分析方法进行对比。