rjags ANOVA with separate variances for each group

Question

I am trying to use rjags to conduct an ANOVA wih separate variances for each of the three groups. I am using the "PlantGrowth" data built into R.

Here is the model:

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))
}

The inits function in this model does not work. Can anyone help?

Answer 1

I would suggest next approach. First, bayesian model needs to be trained. So your steps are fine. Take into account that you are defining extremely low priors on you coefficients. So, when you set your initial values too high, issues related to convergence in the chains can appear. The model chain is fine. I only slightly improved the initial values function. Here the code. In order to complete ANOVA you have to work with all the results for the parameters of your interested. This means computing posterior credible intervals, etc. In that way, your main input will be the parameters from chains. Here the code for obtaining the parameters you want:

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))

The outputs will be stored on d1. Here some rows:

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

With that values you have to work on computing posterior means and posterior credible intervals in order to contrast with the ANOVA method.