1

我正在尝试使用 rjags 对三组中的每一个进行方差分析。我正在使用 R 中内置的“PlantGrowth”数据。

这是模型：

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 函数不起作用。任何人都可以帮忙吗？

4

1 回答 1

1

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

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

使用该值，您必须计算后验均值和后验可信区间，以便与 ANOVA 方法进行对比。

于 2020-09-17T23:04:31.110 回答